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

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

bugfix for commit 3887

  • Property svn:keywords set to Id
File size: 224.9 KB
<
Line 
1!> @file chemistry_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2019 Karlsruhe Institute of Technology
19! Copyright 2017-2019 Freie Universitaet Berlin
20!------------------------------------------------------------------------------!
21!
22! Current revisions:
23! -----------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: chemistry_model_mod.f90 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 SUBROUTINE chem_non_transport_physics_ij( i, j )
2439
2440
2441    INTEGER(iwp), INTENT(IN) ::  i  !< grid index in x-direction
2442    INTEGER(iwp), INTENT(IN) ::  j  !< grid index in y-direction
2443
2444!
2445!-- Calculation of chemical reactions and deposition.
2446    IF ( chem_gasphase_on ) THEN
2447
2448       IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
2449
2450          CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2451          CALL chem_integrate( i, j )
2452          CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2453
2454          IF ( deposition_dry )  THEN
2455             CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2456             CALL chem_depo( i, j )
2457             CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2458          ENDIF
2459
2460       ENDIF
2461
2462    ENDIF
2463
2464 END SUBROUTINE chem_non_transport_physics_ij
2465 
2466!------------------------------------------------------------------------------!
2467! Description:
2468! ------------
2469!> routine for exchange horiz of chemical quantities 
2470!------------------------------------------------------------------------------!
2471 SUBROUTINE chem_exchange_horiz
2472
2473
2474 END SUBROUTINE chem_exchange_horiz
2475
2476 
2477!------------------------------------------------------------------------------!
2478! Description:
2479! ------------
2480!> Subroutine calculating prognostic equations for chemical species
2481!> (vector-optimized).
2482!> Routine is called separately for each chemical species over a loop from
2483!> prognostic_equations.
2484!------------------------------------------------------------------------------!
2485 SUBROUTINE chem_prognostic_equations()
2486
2487
2488    INTEGER ::  i   !< running index
2489    INTEGER ::  j   !< running index
2490    INTEGER ::  k   !< running index
2491
2492    INTEGER(iwp) ::  ilsp   !<
2493
2494
2495    CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2496
2497    DO  ilsp = 1, nspec
2498!
2499!--    Tendency terms for chemical species
2500       tend = 0.0_wp
2501!
2502!--    Advection terms
2503       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2504          IF ( ws_scheme_sca )  THEN
2505             CALL advec_s_ws( chem_species(ilsp)%conc, 'kc' )
2506          ELSE
2507             CALL advec_s_pw( chem_species(ilsp)%conc )
2508          ENDIF
2509       ELSE
2510          CALL advec_s_up( chem_species(ilsp)%conc )
2511       ENDIF
2512!
2513!--    Diffusion terms  (the last three arguments are zero)
2514       CALL diffusion_s( chem_species(ilsp)%conc,                                                  &
2515            surf_def_h(0)%cssws(ilsp,:),                                                           &
2516            surf_def_h(1)%cssws(ilsp,:),                                                           &
2517            surf_def_h(2)%cssws(ilsp,:),                                                           &
2518            surf_lsm_h%cssws(ilsp,:),                                                              &
2519            surf_usm_h%cssws(ilsp,:),                                                              &
2520            surf_def_v(0)%cssws(ilsp,:),                                                           &
2521            surf_def_v(1)%cssws(ilsp,:),                                                           &
2522            surf_def_v(2)%cssws(ilsp,:),                                                           &
2523            surf_def_v(3)%cssws(ilsp,:),                                                           &
2524            surf_lsm_v(0)%cssws(ilsp,:),                                                           &
2525            surf_lsm_v(1)%cssws(ilsp,:),                                                           &
2526            surf_lsm_v(2)%cssws(ilsp,:),                                                           &
2527            surf_lsm_v(3)%cssws(ilsp,:),                                                           &
2528            surf_usm_v(0)%cssws(ilsp,:),                                                           &
2529            surf_usm_v(1)%cssws(ilsp,:),                                                           &
2530            surf_usm_v(2)%cssws(ilsp,:),                                                           &
2531            surf_usm_v(3)%cssws(ilsp,:) )
2532!
2533!--    Prognostic equation for chemical species
2534       DO  i = nxl, nxr
2535          DO  j = nys, nyn
2536             DO  k = nzb+1, nzt
2537                chem_species(ilsp)%conc_p(k,j,i) =   chem_species(ilsp)%conc(k,j,i)                &
2538                     + ( dt_3d  *                                                                  &
2539                     (   tsc(2) * tend(k,j,i)                                                      &
2540                     + tsc(3) * chem_species(ilsp)%tconc_m(k,j,i)                                  &
2541                     )                                                                             &
2542                     - tsc(5) * rdf_sc(k)                                                          &
2543                     * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )     &
2544                     )                                                                             &
2545                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2546
2547                IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
2548                   chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)
2549                ENDIF
2550             ENDDO
2551          ENDDO
2552       ENDDO
2553!
2554!--    Calculate tendencies for the next Runge-Kutta step
2555       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2556          IF ( intermediate_timestep_count == 1 )  THEN
2557             DO  i = nxl, nxr
2558                DO  j = nys, nyn
2559                   DO  k = nzb+1, nzt
2560                      chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
2561                   ENDDO
2562                ENDDO
2563             ENDDO
2564          ELSEIF ( intermediate_timestep_count < &
2565               intermediate_timestep_count_max )  THEN
2566             DO  i = nxl, nxr
2567                DO  j = nys, nyn
2568                   DO  k = nzb+1, nzt
2569                      chem_species(ilsp)%tconc_m(k,j,i) = - 9.5625_wp * tend(k,j,i)                &
2570                           + 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
2571                   ENDDO
2572                ENDDO
2573             ENDDO
2574          ENDIF
2575       ENDIF
2576
2577    ENDDO
2578
2579    CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )
2580
2581 END SUBROUTINE chem_prognostic_equations
2582
2583
2584!------------------------------------------------------------------------------!
2585! Description:
2586! ------------
2587!> Subroutine calculating prognostic equations for chemical species
2588!> (cache-optimized).
2589!> Routine is called separately for each chemical species over a loop from
2590!> prognostic_equations.
2591!------------------------------------------------------------------------------!
2592 SUBROUTINE chem_prognostic_equations_ij( i, j, i_omp_start, tn )
2593
2594
2595    INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn
2596    INTEGER(iwp) :: ilsp
2597!
2598!-- local variables
2599
2600    INTEGER :: k
2601
2602    DO  ilsp = 1, nspec
2603!
2604!--    Tendency-terms for chem spcs.
2605       tend(:,j,i) = 0.0_wp
2606!
2607!--    Advection terms
2608       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2609          IF ( ws_scheme_sca )  THEN
2610             CALL advec_s_ws( i, j, chem_species(ilsp)%conc, 'kc', chem_species(ilsp)%flux_s_cs,   &
2611                              chem_species(ilsp)%diss_s_cs, chem_species(ilsp)%flux_l_cs,          &
2612                              chem_species(ilsp)%diss_l_cs, i_omp_start, tn )
2613          ELSE
2614             CALL advec_s_pw( i, j, chem_species(ilsp)%conc )
2615          ENDIF
2616       ELSE
2617          CALL advec_s_up( i, j, chem_species(ilsp)%conc )
2618       ENDIF
2619!
2620!--    Diffusion terms (the last three arguments are zero)
2621
2622       CALL diffusion_s( i, j, chem_species(ilsp)%conc,                                            &
2623            surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),                              &
2624            surf_def_h(2)%cssws(ilsp,:),                                                           &
2625            surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:),                                    &
2626            surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),                              &
2627            surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),                              &
2628            surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:),                              &
2629            surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:),                              &
2630            surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),                              &
2631            surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
2632!
2633!--    Prognostic equation for chem spcs
2634       DO  k = nzb+1, nzt
2635          chem_species(ilsp)%conc_p(k,j,i) = chem_species(ilsp)%conc(k,j,i) + ( dt_3d  *           &
2636               ( tsc(2) * tend(k,j,i) +                                                            &
2637               tsc(3) * chem_species(ilsp)%tconc_m(k,j,i) )                                        &
2638               - tsc(5) * rdf_sc(k)                                                                &
2639               * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )           &
2640               )                                                                                   &
2641               * MERGE( 1.0_wp, 0.0_wp,                                                            &
2642               BTEST( wall_flags_0(k,j,i), 0 )                                                     &
2643               )
2644
2645          IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
2646             chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6
2647          ENDIF
2648       ENDDO
2649!
2650!--    Calculate tendencies for the next Runge-Kutta step
2651       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2652          IF ( intermediate_timestep_count == 1 )  THEN
2653             DO  k = nzb+1, nzt
2654                chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
2655             ENDDO
2656          ELSEIF ( intermediate_timestep_count < &
2657               intermediate_timestep_count_max )  THEN
2658             DO  k = nzb+1, nzt
2659                chem_species(ilsp)%tconc_m(k,j,i) = -9.5625_wp * tend(k,j,i) +                     &
2660                     5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
2661             ENDDO
2662          ENDIF
2663       ENDIF
2664
2665    ENDDO
2666
2667 END SUBROUTINE chem_prognostic_equations_ij
2668
2669
2670!------------------------------------------------------------------------------!
2671! Description:
2672! ------------
2673!> Subroutine to read restart data of chemical species
2674!------------------------------------------------------------------------------!
2675 SUBROUTINE chem_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
2676                            nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,   &
2677                            nys_on_file, tmp_3d, found )
2678
2679    USE control_parameters
2680
2681
2682    CHARACTER (LEN=20) :: spc_name_av !<   
2683
2684    INTEGER(iwp) ::  lsp             !<
2685    INTEGER(iwp) ::  k               !<
2686    INTEGER(iwp) ::  nxlc            !<
2687    INTEGER(iwp) ::  nxlf            !<
2688    INTEGER(iwp) ::  nxl_on_file     !<   
2689    INTEGER(iwp) ::  nxrc            !<
2690    INTEGER(iwp) ::  nxrf            !<
2691    INTEGER(iwp) ::  nxr_on_file     !<   
2692    INTEGER(iwp) ::  nync            !<
2693    INTEGER(iwp) ::  nynf            !<
2694    INTEGER(iwp) ::  nyn_on_file     !<   
2695    INTEGER(iwp) ::  nysc            !<
2696    INTEGER(iwp) ::  nysf            !<
2697    INTEGER(iwp) ::  nys_on_file     !<   
2698
2699    LOGICAL, INTENT(OUT) :: found
2700
2701    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
2702
2703
2704    found = .FALSE. 
2705
2706
2707    IF ( ALLOCATED(chem_species) )  THEN
2708
2709       DO  lsp = 1, nspec
2710
2711          !< for time-averaged chemical conc.
2712          spc_name_av  =  TRIM( chem_species(lsp)%name )//'_av'
2713
2714          IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) )    &
2715             THEN
2716             !< read data into tmp_3d
2717             IF ( k == 1 )  READ ( 13 )  tmp_3d 
2718             !< fill ..%conc in the restart run   
2719             chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                    &
2720                  nxlc-nbgp:nxrc+nbgp) =                  & 
2721                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2722             found = .TRUE.
2723          ELSEIF (restart_string(1:length) == spc_name_av )  THEN
2724             IF ( k == 1 )  READ ( 13 )  tmp_3d
2725             chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,                 &
2726                  nxlc-nbgp:nxrc+nbgp) =               &
2727                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2728             found = .TRUE.
2729          ENDIF
2730
2731       ENDDO
2732
2733    ENDIF
2734
2735
2736 END SUBROUTINE chem_rrd_local
2737
2738
2739!-------------------------------------------------------------------------------!
2740!> Description:
2741!> Calculation of horizontally averaged profiles
2742!> This routine is called for every statistic region (sr) defined by the user,
2743!> but at least for the region "total domain" (sr=0).
2744!> quantities.
2745!-------------------------------------------------------------------------------!
2746 SUBROUTINE chem_statistics( mode, sr, tn )
2747
2748
2749    USE arrays_3d
2750
2751    USE statistics
2752
2753
2754    CHARACTER (LEN=*) ::  mode   !<
2755
2756    INTEGER(iwp) ::  i    !< running index on x-axis
2757    INTEGER(iwp) ::  j    !< running index on y-axis
2758    INTEGER(iwp) ::  k    !< vertical index counter
2759    INTEGER(iwp) ::  sr   !< statistical region
2760    INTEGER(iwp) ::  tn   !< thread number
2761    INTEGER(iwp) ::  lpr  !< running index chem spcs
2762
2763    IF ( mode == 'profiles' )  THEN
2764       !
2765!
2766!--    Sample on how to calculate horizontally averaged profiles of user-
2767!--    defined quantities. Each quantity is identified by the index
2768!--    "pr_palm+#" where "#" is an integer starting from 1. These
2769!--    user-profile-numbers must also be assigned to the respective strings
2770!--    given by data_output_pr_cs in routine user_check_data_output_pr.
2771!--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
2772!--                     w*pt*), dim-4 = statistical region.
2773
2774!$OMP DO
2775       DO  i = nxl, nxr
2776          DO  j = nys, nyn
2777             DO  k = nzb, nzt+1
2778                DO lpr = 1, cs_pr_count
2779
2780                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
2781                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
2782                        rmask(j,i,sr)  *                                   &
2783                        MERGE( 1.0_wp, 0.0_wp,                             &
2784                        BTEST( wall_flags_0(k,j,i), 22 ) )
2785                ENDDO
2786             ENDDO
2787          ENDDO
2788       ENDDO
2789    ELSEIF ( mode == 'time_series' )  THEN
2790!      @todo
2791    ENDIF
2792
2793 END SUBROUTINE chem_statistics
2794
2795
2796!------------------------------------------------------------------------------!
2797! Description:
2798! ------------
2799!> Subroutine for swapping of timelevels for chemical species
2800!> called out from subroutine swap_timelevel
2801!------------------------------------------------------------------------------!
2802
2803
2804 SUBROUTINE chem_swap_timelevel( level )
2805
2806
2807    INTEGER(iwp), INTENT(IN) ::  level
2808!
2809!-- local variables
2810    INTEGER(iwp)             ::  lsp
2811
2812
2813    IF ( level == 0 )  THEN
2814       DO  lsp=1, nvar                                       
2815          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
2816          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
2817       ENDDO
2818    ELSE
2819       DO  lsp=1, nvar                                       
2820          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
2821          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
2822       ENDDO
2823    ENDIF
2824
2825    RETURN
2826 END SUBROUTINE chem_swap_timelevel
2827
2828
2829!------------------------------------------------------------------------------!
2830! Description:
2831! ------------
2832!> Subroutine to write restart data for chemistry model
2833!------------------------------------------------------------------------------!
2834 SUBROUTINE chem_wrd_local
2835
2836
2837    INTEGER(iwp) ::  lsp  !< running index for chem spcs.
2838
2839    DO  lsp = 1, nspec
2840       CALL wrd_write_string( TRIM( chem_species(lsp)%name ) )
2841       WRITE ( 14 )  chem_species(lsp)%conc
2842       CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
2843       WRITE ( 14 )  chem_species(lsp)%conc_av
2844    ENDDO
2845
2846 END SUBROUTINE chem_wrd_local
2847
2848
2849!!! sB remove blanks again later !!!
2850
2851
2852
2853
2854
2855
2856
2857
2858
2859
2860
2861
2862
2863
2864
2865
2866
2867
2868
2869
2870
2871
2872
2873
2874
2875
2876
2877
2878
2879
2880
2881
2882
2883
2884
2885
2886
2887
2888!-------------------------------------------------------------------------------!
2889! Description:
2890! ------------
2891!> Subroutine to calculate the deposition of gases and PMs. For now deposition
2892!> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
2893!> considered. The deposition of particles is derived following Zhang et al.,
2894!> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
2895!>     
2896!> @TODO: Consider deposition on vertical surfaces   
2897!> @TODO: Consider overlaying horizontal surfaces
2898!> @TODO: Consider resolved vegetation   
2899!> @TODO: Check error messages
2900!-------------------------------------------------------------------------------!
2901 SUBROUTINE chem_depo( i, j )
2902
2903    USE control_parameters,                                                 &   
2904         ONLY:  dt_3d, intermediate_timestep_count, latitude
2905
2906    USE arrays_3d,                                                          &
2907         ONLY:  dzw, rho_air_zw
2908
2909    USE date_and_time_mod,                                                  &
2910         ONLY:  day_of_year
2911
2912    USE surface_mod,                                                        &
2913         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h,        &
2914         surf_type, surf_usm_h
2915
2916    USE radiation_model_mod,                                                &
2917         ONLY:  cos_zenith
2918
2919
2920    INTEGER(iwp), INTENT(IN) ::  i
2921    INTEGER(iwp), INTENT(IN) ::  j
2922    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
2923    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
2924    INTEGER(iwp) ::  lu_palm                       !< index of PALM LSM vegetation_type at current surface element
2925    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
2926    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current surface element
2927    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current surface element
2928    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current surface element
2929    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface element
2930    INTEGER(iwp) ::  lu_dep                        !< matching DEPAC LU to lu_palm
2931    INTEGER(iwp) ::  lup_dep                       !< matching DEPAC LU to lup_palm
2932    INTEGER(iwp) ::  luw_dep                       !< matching DEPAC LU to luw_palm
2933    INTEGER(iwp) ::  luu_dep                       !< matching DEPAC LU to luu_palm
2934    INTEGER(iwp) ::  lug_dep                       !< matching DEPAC LU to lug_palm
2935    INTEGER(iwp) ::  lud_dep                       !< matching DEPAC LU to lud_palm
2936    INTEGER(iwp) ::  m                             !< index for horizontal surfaces
2937
2938    INTEGER(iwp) ::  pspec                         !< running index
2939    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
2940!
2941!-- Vegetation                                               !< Assign PALM classes to DEPAC land use classes
2942    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM
2943    INTEGER(iwp) ::  ind_luv_b_soil = 1                      !<  assigned to ilu_desert
2944    INTEGER(iwp) ::  ind_luv_mixed_crops = 2                 !<  assigned to ilu_arable
2945    INTEGER(iwp) ::  ind_luv_s_grass = 3                     !<  assigned to ilu_grass
2946    INTEGER(iwp) ::  ind_luv_ev_needle_trees = 4             !<  assigned to ilu_coniferous_forest
2947    INTEGER(iwp) ::  ind_luv_de_needle_trees = 5             !<  assigned to ilu_coniferous_forest
2948    INTEGER(iwp) ::  ind_luv_ev_broad_trees = 6              !<  assigned to ilu_tropical_forest
2949    INTEGER(iwp) ::  ind_luv_de_broad_trees = 7              !<  assigned to ilu_deciduous_forest
2950    INTEGER(iwp) ::  ind_luv_t_grass = 8                     !<  assigned to ilu_grass
2951    INTEGER(iwp) ::  ind_luv_desert = 9                      !<  assigned to ilu_desert
2952    INTEGER(iwp) ::  ind_luv_tundra = 10                     !<  assigned to ilu_other
2953    INTEGER(iwp) ::  ind_luv_irr_crops = 11                  !<  assigned to ilu_arable
2954    INTEGER(iwp) ::  ind_luv_semidesert = 12                 !<  assigned to ilu_other
2955    INTEGER(iwp) ::  ind_luv_ice = 13                        !<  assigned to ilu_ice
2956    INTEGER(iwp) ::  ind_luv_marsh = 14                      !<  assigned to ilu_other
2957    INTEGER(iwp) ::  ind_luv_ev_shrubs = 15                  !<  assigned to ilu_mediterrean_scrub
2958    INTEGER(iwp) ::  ind_luv_de_shrubs = 16                  !<  assigned to ilu_mediterrean_scrub
2959    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest (ave(decid+conif))
2960    INTEGER(iwp) ::  ind_luv_intrup_forest = 18              !<  assigned to ilu_other (ave(other+decid))
2961!
2962!-- Water
2963    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM 
2964    INTEGER(iwp) ::  ind_luw_lake = 1                        !<  assigned to ilu_water_inland
2965    INTEGER(iwp) ::  ind_luw_river = 2                       !<  assigned to ilu_water_inland
2966    INTEGER(iwp) ::  ind_luw_ocean = 3                       !<  assigned to ilu_water_sea
2967    INTEGER(iwp) ::  ind_luw_pond = 4                        !<  assigned to ilu_water_inland
2968    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland
2969!
2970!-- Pavement
2971    INTEGER(iwp) ::  ind_lup_user = 0                        !<  ERROR as no class given in PALM
2972    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !<  assigned to ilu_desert
2973    INTEGER(iwp) ::  ind_lup_asph = 2                        !<  assigned to ilu_desert
2974    INTEGER(iwp) ::  ind_lup_conc = 3                        !<  assigned to ilu_desert
2975    INTEGER(iwp) ::  ind_lup_sett = 4                        !<  assigned to ilu_desert
2976    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert
2977    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !<  assigned to ilu_desert
2978    INTEGER(iwp) ::  ind_lup_metal = 7                       !<  assigned to ilu_desert
2979    INTEGER(iwp) ::  ind_lup_wood = 8                        !<  assigned to ilu_desert
2980    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert
2981    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !<  assigned to ilu_desert
2982    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !<  assigned to ilu_desert
2983    INTEGER(iwp) ::  ind_lup_woodchips = 12                  !<  assigned to ilu_desert
2984    INTEGER(iwp) ::  ind_lup_tartan = 13                     !<  assigned to ilu_desert
2985    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !<  assigned to ilu_desert
2986    INTEGER(iwp) ::  ind_lup_clay = 15                       !<  assigned to ilu_desert
2987!
2988!-- Particle parameters according to the respective aerosol classes (PM25, PM10)
2989    INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
2990    INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
2991    INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
2992
2993    INTEGER(iwp) ::  part_type          !< index for particle type (PM10 or PM25) in particle_pars
2994
2995    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
2996
2997    REAL(wp) ::  dt_chem                !< length of chem time step
2998    REAL(wp) ::  dh                     !< vertical grid size
2999    REAL(wp) ::  inv_dh                 !< inverse of vertical grid size
3000    REAL(wp) ::  dt_dh                  !< dt_chem/dh
3001
3002    REAL(wp) ::  dens              !< density at layer k at i,j 
3003    REAL(wp) ::  r_aero_surf       !< aerodynamic resistance (s/m) at current surface element
3004    REAL(wp) ::  ustar_surf        !< ustar at current surface element
3005    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element
3006    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface element
3007    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
3008    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
3009    REAL(wp) ::  lai               !< leaf area index at current surface element
3010    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
3011
3012    REAL(wp) ::  slinnfac       
3013    REAL(wp) ::  visc              !< Viscosity
3014    REAL(wp) ::  vs                !< Sedimentation velocity
3015    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
3016    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
3017    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
3018    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
3019
3020    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
3021    REAL(wp) ::  diffusivity       !< diffusivity
3022
3023
3024    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface element
3025    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface element
3026    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface element
3027    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface element
3028    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface element
3029    REAL(wp), DIMENSION(nspec) ::  bud_lud      !< budget for USM windows at current surface element
3030    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
3031    REAL(wp), DIMENSION(nspec) ::  conc_ijk     !< concentration at i,j,k
3032    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to zero for all species!
3033                                                 
3034
3035    REAL(wp) ::  temp_tmp       !< temperatur at i,j,k
3036    REAL(wp) ::  ts             !< surface temperatur in degrees celsius
3037    REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
3038!
3039!-- Particle parameters (PM10 (1), PM25 (2))
3040!-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
3041!-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
3042    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/ &
3043         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &  !<  1
3044         0.7e-6_wp, 1.14e3_wp, 1.082_wp &   !<  2
3045         /), (/ 3, 2 /) )
3046
3047    LOGICAL ::  match_lsm     !< flag indicating natural-type surface
3048    LOGICAL ::  match_usm     !< flag indicating urban-type surface
3049!
3050!-- List of names of possible tracers
3051    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/ &
3052         'NO2           ', &    !< NO2
3053         'NO            ', &    !< NO
3054         'O3            ', &    !< O3
3055         'CO            ', &    !< CO
3056         'form          ', &    !< FORM
3057         'ald           ', &    !< ALD
3058         'pan           ', &    !< PAN
3059         'mgly          ', &    !< MGLY
3060         'par           ', &    !< PAR
3061         'ole           ', &    !< OLE
3062         'eth           ', &    !< ETH
3063         'tol           ', &    !< TOL
3064         'cres          ', &    !< CRES
3065         'xyl           ', &    !< XYL
3066         'SO4a_f        ', &    !< SO4a_f
3067         'SO2           ', &    !< SO2
3068         'HNO2          ', &    !< HNO2
3069         'CH4           ', &    !< CH4
3070         'NH3           ', &    !< NH3
3071         'NO3           ', &    !< NO3
3072         'OH            ', &    !< OH
3073         'HO2           ', &    !< HO2
3074         'N2O5          ', &    !< N2O5
3075         'SO4a_c        ', &    !< SO4a_c
3076         'NH4a_f        ', &    !< NH4a_f
3077         'NO3a_f        ', &    !< NO3a_f
3078         'NO3a_c        ', &    !< NO3a_c
3079         'C2O3          ', &    !< C2O3
3080         'XO2           ', &    !< XO2
3081         'XO2N          ', &    !< XO2N
3082         'cro           ', &    !< CRO
3083         'HNO3          ', &    !< HNO3
3084         'H2O2          ', &    !< H2O2
3085         'iso           ', &    !< ISO
3086         'ispd          ', &    !< ISPD
3087         'to2           ', &    !< TO2
3088         'open          ', &    !< OPEN
3089         'terp          ', &    !< TERP
3090         'ec_f          ', &    !< EC_f
3091         'ec_c          ', &    !< EC_c
3092         'pom_f         ', &    !< POM_f
3093         'pom_c         ', &    !< POM_c
3094         'ppm_f         ', &    !< PPM_f
3095         'ppm_c         ', &    !< PPM_c
3096         'na_ff         ', &    !< Na_ff
3097         'na_f          ', &    !< Na_f
3098         'na_c          ', &    !< Na_c
3099         'na_cc         ', &    !< Na_cc
3100         'na_ccc        ', &    !< Na_ccc
3101         'dust_ff       ', &    !< dust_ff
3102         'dust_f        ', &    !< dust_f
3103         'dust_c        ', &    !< dust_c
3104         'dust_cc       ', &    !< dust_cc
3105         'dust_ccc      ', &    !< dust_ccc
3106         'tpm10         ', &    !< tpm10
3107         'tpm25         ', &    !< tpm25
3108         'tss           ', &    !< tss
3109         'tdust         ', &    !< tdust
3110         'tc            ', &    !< tc
3111         'tcg           ', &    !< tcg
3112         'tsoa          ', &    !< tsoa
3113         'tnmvoc        ', &    !< tnmvoc
3114         'SOxa          ', &    !< SOxa
3115         'NOya          ', &    !< NOya
3116         'NHxa          ', &    !< NHxa
3117         'NO2_obs       ', &    !< NO2_obs
3118         'tpm10_biascorr', &    !< tpm10_biascorr
3119         'tpm25_biascorr', &    !< tpm25_biascorr
3120         'O3_biascorr   ' /)    !< o3_biascorr
3121!
3122!-- tracer mole mass:
3123    REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
3124         xm_O * 2 + xm_N, &                         !< NO2
3125         xm_O + xm_N, &                             !< NO
3126         xm_O * 3, &                                !< O3
3127         xm_C + xm_O, &                             !< CO
3128         xm_H * 2 + xm_C + xm_O, &                  !< FORM
3129         xm_H * 3 + xm_C * 2 + xm_O, &              !< ALD
3130         xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, &   !< PAN
3131         xm_H * 4 + xm_C * 3 + xm_O * 2, &          !< MGLY
3132         xm_H * 3 + xm_C, &                         !< PAR
3133         xm_H * 3 + xm_C * 2, &                     !< OLE
3134         xm_H * 4 + xm_C * 2, &                     !< ETH
3135         xm_H * 8 + xm_C * 7, &                     !< TOL
3136         xm_H * 8 + xm_C * 7 + xm_O, &              !< CRES
3137         xm_H * 10 + xm_C * 8, &                    !< XYL
3138         xm_S + xm_O * 4, &                         !< SO4a_f
3139         xm_S + xm_O * 2, &                         !< SO2
3140         xm_H + xm_O * 2 + xm_N, &                  !< HNO2
3141         xm_H * 4 + xm_C, &                         !< CH4
3142         xm_H * 3 + xm_N, &                         !< NH3
3143         xm_O * 3 + xm_N, &                         !< NO3
3144         xm_H + xm_O, &                             !< OH
3145         xm_H + xm_O * 2, &                         !< HO2
3146         xm_O * 5 + xm_N * 2, &                     !< N2O5
3147         xm_S + xm_O * 4, &                         !< SO4a_c
3148         xm_H * 4 + xm_N, &                         !< NH4a_f
3149         xm_O * 3 + xm_N, &                         !< NO3a_f
3150         xm_O * 3 + xm_N, &                         !< NO3a_c
3151         xm_C * 2 + xm_O * 3, &                     !< C2O3
3152         xm_dummy, &                                !< XO2
3153         xm_dummy, &                                !< XO2N
3154         xm_dummy, &                                !< CRO
3155         xm_H + xm_O * 3 + xm_N, &                  !< HNO3
3156         xm_H * 2 + xm_O * 2, &                     !< H2O2
3157         xm_H * 8 + xm_C * 5, &                     !< ISO
3158         xm_dummy, &                                !< ISPD
3159         xm_dummy, &                                !< TO2
3160         xm_dummy, &                                !< OPEN
3161         xm_H * 16 + xm_C * 10, &                   !< TERP
3162         xm_dummy, &                                !< EC_f
3163         xm_dummy, &                                !< EC_c
3164         xm_dummy, &                                !< POM_f
3165         xm_dummy, &                                !< POM_c
3166         xm_dummy, &                                !< PPM_f
3167         xm_dummy, &                                !< PPM_c
3168         xm_Na, &                                   !< Na_ff
3169         xm_Na, &                                   !< Na_f
3170         xm_Na, &                                   !< Na_c
3171         xm_Na, &                                   !< Na_cc
3172         xm_Na, &                                   !< Na_ccc
3173         xm_dummy, &                                !< dust_ff
3174         xm_dummy, &                                !< dust_f
3175         xm_dummy, &                                !< dust_c
3176         xm_dummy, &                                !< dust_cc
3177         xm_dummy, &                                !< dust_ccc
3178         xm_dummy, &                                !< tpm10
3179         xm_dummy, &                                !< tpm25
3180         xm_dummy, &                                !< tss
3181         xm_dummy, &                                !< tdust
3182         xm_dummy, &                                !< tc
3183         xm_dummy, &                                !< tcg
3184         xm_dummy, &                                !< tsoa
3185         xm_dummy, &                                !< tnmvoc
3186         xm_dummy, &                                !< SOxa
3187         xm_dummy, &                                !< NOya
3188         xm_dummy, &                                !< NHxa
3189         xm_O * 2 + xm_N, &                         !< NO2_obs
3190         xm_dummy, &                                !< tpm10_biascorr
3191         xm_dummy, &                                !< tpm25_biascorr
3192         xm_O * 3 /)                                !< o3_biascorr
3193!
3194!-- Initialize surface element m
3195    m = 0
3196    k = 0
3197!
3198!-- LSM or USM surface present at i,j:
3199!-- Default surfaces are NOT considered for deposition
3200    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
3201    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
3202!
3203!--For LSM surfaces
3204
3205    IF ( match_lsm )  THEN
3206!
3207!--    Get surface element information at i,j:
3208       m = surf_lsm_h%start_index(j,i)
3209       k = surf_lsm_h%k(m)
3210!
3211!--    Get needed variables for surface element m
3212       ustar_surf  = surf_lsm_h%us(m)
3213       z0h_surf    = surf_lsm_h%z0h(m)
3214       r_aero_surf = surf_lsm_h%r_a(m)
3215       solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
3216       lai = surf_lsm_h%lai(m)
3217       sai = lai + 1
3218!
3219!--    For small grid spacing neglect R_a
3220       IF ( dzw(k) <= 1.0 )  THEN
3221          r_aero_surf = 0.0_wp
3222       ENDIF
3223!
3224!--    Initialize lu's
3225       lu_palm = 0
3226       lu_dep = 0
3227       lup_palm = 0
3228       lup_dep = 0
3229       luw_palm = 0
3230       luw_dep = 0
3231!
3232!--    Initialize budgets
3233       bud_luv = 0.0_wp
3234       bud_lup = 0.0_wp
3235       bud_luw = 0.0_wp
3236!
3237!--    Get land use for i,j and assign to DEPAC lu
3238       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
3239          lu_palm = surf_lsm_h%vegetation_type(m)
3240          IF ( lu_palm == ind_luv_user )  THEN
3241             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
3242             CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
3243          ELSEIF ( lu_palm == ind_luv_b_soil )  THEN
3244             lu_dep = 9
3245          ELSEIF ( lu_palm == ind_luv_mixed_crops )  THEN
3246             lu_dep = 2
3247          ELSEIF ( lu_palm == ind_luv_s_grass )  THEN
3248             lu_dep = 1
3249          ELSEIF ( lu_palm == ind_luv_ev_needle_trees )  THEN
3250             lu_dep = 4
3251          ELSEIF ( lu_palm == ind_luv_de_needle_trees )  THEN
3252             lu_dep = 4
3253          ELSEIF ( lu_palm == ind_luv_ev_broad_trees )  THEN
3254             lu_dep = 12
3255          ELSEIF ( lu_palm == ind_luv_de_broad_trees )  THEN
3256             lu_dep = 5
3257          ELSEIF ( lu_palm == ind_luv_t_grass )  THEN
3258             lu_dep = 1
3259          ELSEIF ( lu_palm == ind_luv_desert )  THEN
3260             lu_dep = 9
3261          ELSEIF ( lu_palm == ind_luv_tundra )  THEN
3262             lu_dep = 8
3263          ELSEIF ( lu_palm == ind_luv_irr_crops )  THEN
3264             lu_dep = 2
3265          ELSEIF ( lu_palm == ind_luv_semidesert )  THEN
3266             lu_dep = 8
3267          ELSEIF ( lu_palm == ind_luv_ice )  THEN
3268             lu_dep = 10
3269          ELSEIF ( lu_palm == ind_luv_marsh )  THEN
3270             lu_dep = 8
3271          ELSEIF ( lu_palm == ind_luv_ev_shrubs )  THEN
3272             lu_dep = 14
3273          ELSEIF ( lu_palm == ind_luv_de_shrubs )  THEN
3274             lu_dep = 14
3275          ELSEIF ( lu_palm == ind_luv_mixed_forest )  THEN
3276             lu_dep = 4
3277          ELSEIF ( lu_palm == ind_luv_intrup_forest )  THEN
3278             lu_dep = 8     
3279          ENDIF
3280       ENDIF
3281
3282       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
3283          lup_palm = surf_lsm_h%pavement_type(m)
3284          IF ( lup_palm == ind_lup_user )  THEN
3285             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3286             CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 )
3287          ELSEIF ( lup_palm == ind_lup_asph_conc )  THEN
3288             lup_dep = 9
3289          ELSEIF ( lup_palm == ind_lup_asph )  THEN
3290             lup_dep = 9
3291          ELSEIF ( lup_palm ==  ind_lup_conc )  THEN
3292             lup_dep = 9
3293          ELSEIF ( lup_palm ==  ind_lup_sett )  THEN
3294             lup_dep = 9
3295          ELSEIF ( lup_palm == ind_lup_pav_stones )  THEN
3296             lup_dep = 9
3297          ELSEIF ( lup_palm == ind_lup_cobblest )  THEN
3298             lup_dep = 9       
3299          ELSEIF ( lup_palm == ind_lup_metal )  THEN
3300             lup_dep = 9
3301          ELSEIF ( lup_palm == ind_lup_wood )  THEN
3302             lup_dep = 9   
3303          ELSEIF ( lup_palm == ind_lup_gravel )  THEN
3304             lup_dep = 9
3305          ELSEIF ( lup_palm == ind_lup_f_gravel )  THEN
3306             lup_dep = 9
3307          ELSEIF ( lup_palm == ind_lup_pebblest )  THEN
3308             lup_dep = 9
3309          ELSEIF ( lup_palm == ind_lup_woodchips )  THEN
3310             lup_dep = 9
3311          ELSEIF ( lup_palm == ind_lup_tartan )  THEN
3312             lup_dep = 9
3313          ELSEIF ( lup_palm == ind_lup_art_turf )  THEN
3314             lup_dep = 9
3315          ELSEIF ( lup_palm == ind_lup_clay )  THEN
3316             lup_dep = 9
3317          ENDIF
3318       ENDIF
3319
3320       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
3321          luw_palm = surf_lsm_h%water_type(m)     
3322          IF ( luw_palm == ind_luw_user )  THEN
3323             message_string = 'No water type defined. Please define water type to enable deposition calculation'
3324             CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 )
3325          ELSEIF ( luw_palm ==  ind_luw_lake )  THEN
3326             luw_dep = 13
3327          ELSEIF ( luw_palm == ind_luw_river )  THEN
3328             luw_dep = 13
3329          ELSEIF ( luw_palm == ind_luw_ocean )  THEN
3330             luw_dep = 6
3331          ELSEIF ( luw_palm == ind_luw_pond )  THEN
3332             luw_dep = 13
3333          ELSEIF ( luw_palm == ind_luw_fountain )  THEN
3334             luw_dep = 13 
3335          ENDIF
3336       ENDIF
3337!
3338!--    Set wetness indicator to dry or wet for lsm vegetation or pavement
3339       IF ( surf_lsm_h%c_liq(m) > 0 )  THEN
3340          nwet = 1
3341       ELSE
3342          nwet = 0
3343       ENDIF
3344!
3345!--    Compute length of time step
3346       IF ( call_chem_at_all_substeps )  THEN
3347          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
3348       ELSE
3349          dt_chem = dt_3d
3350       ENDIF
3351
3352       dh = dzw(k)
3353       inv_dh = 1.0_wp / dh
3354       dt_dh = dt_chem/dh
3355!
3356!--    Concentration at i,j,k
3357       DO  lsp = 1, nspec
3358          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
3359       ENDDO
3360
3361!--    Temperature at i,j,k
3362       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
3363
3364       ts       = temp_tmp - 273.15  !< in degrees celcius
3365!
3366!--    Viscosity of air
3367       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
3368!
3369!--    Air density at k
3370       dens = rho_air_zw(k)
3371!
3372!--    Calculate relative humidity from specific humidity for DEPAC
3373       qv_tmp = MAX(q(k,j,i),0.0_wp)
3374       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
3375!
3376!-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
3377!-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
3378!
3379!--    Vegetation
3380       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
3381
3382          slinnfac = 1.0_wp
3383!
3384!--       Get deposition velocity vd
3385          DO  lsp = 1, nvar
3386!
3387!--          Initialize
3388             vs = 0.0_wp
3389             vd_lu = 0.0_wp
3390             rs = 0.0_wp
3391             rb = 0.0_wp
3392             rc_tot = 0.0_wp
3393             IF ( spc_names(lsp) == 'PM10' )  THEN
3394                part_type = 1
3395!
3396!--             Sedimentation velocity
3397                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3398                     particle_pars(ind_p_size, part_type), &
3399                     particle_pars(ind_p_slip, part_type), &
3400                     visc)
3401
3402                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3403                     vs, &
3404                     particle_pars(ind_p_size, part_type), &
3405                     particle_pars(ind_p_slip, part_type), &
3406                     nwet, temp_tmp, dens, visc, &
3407                     lu_dep,  &
3408                     r_aero_surf, ustar_surf )
3409
3410                bud_luv(lsp) = - conc_ijk(lsp) * &
3411                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3412
3413
3414             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3415                part_type = 2
3416!
3417!--             Sedimentation velocity
3418                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3419                     particle_pars(ind_p_size, part_type), &
3420                     particle_pars(ind_p_slip, part_type), &
3421                     visc)
3422
3423                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3424                     vs, &
3425                     particle_pars(ind_p_size, part_type), &
3426                     particle_pars(ind_p_slip, part_type), &
3427                     nwet, temp_tmp, dens, visc, &
3428                     lu_dep , &
3429                     r_aero_surf, ustar_surf )
3430
3431                bud_luv(lsp) = - conc_ijk(lsp) * &
3432                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3433
3434             ELSE !< GASES
3435!
3436!--             Read spc_name of current species for gas parameter
3437                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
3438                   i_pspec = 0
3439                   DO  pspec = 1, nposp
3440                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3441                         i_pspec = pspec
3442                      END IF
3443                   ENDDO
3444
3445                ELSE
3446!
3447!--             For now species not deposited
3448                   CYCLE
3449                ENDIF
3450!
3451!--             Factor used for conversion from ppb to ug/m3 :
3452!--             ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3453!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3454!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
3455!--             thus:
3456!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3457!--             Use density at k:
3458
3459                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3460!
3461!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3462                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3463                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
3464!
3465!--             Diffusivity for DEPAC relevant gases
3466!--             Use default value
3467                diffusivity            = 0.11e-4
3468!
3469!--             overwrite with known coefficients of diffusivity from Massman (1998)
3470                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3471                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3472                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3473                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3474                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3475                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3476                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3477!
3478!--             Get quasi-laminar boundary layer resistance rb:
3479                CALL get_rb_cell( (lu_dep == ilu_water_sea) .OR. (lu_dep == ilu_water_inland), &
3480                     z0h_surf, ustar_surf, diffusivity, &
3481                     rb )
3482!
3483!--             Get rc_tot
3484                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, solar_rad, cos_zenith, &
3485                     rh_surf, lai, sai, nwet, lu_dep, 2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
3486                     r_aero_surf , rb )
3487!
3488!--             Calculate budget
3489                IF ( rc_tot <= 0.0 )  THEN
3490
3491                   bud_luv(lsp) = 0.0_wp
3492
3493                ELSE
3494
3495                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3496
3497                   bud_luv(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3498                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3499                ENDIF
3500
3501             ENDIF
3502          ENDDO
3503       ENDIF
3504!
3505!--    Pavement
3506       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
3507!
3508!--       No vegetation on pavements:
3509          lai = 0.0_wp
3510          sai = 0.0_wp
3511
3512          slinnfac = 1.0_wp
3513!
3514!--       Get vd
3515          DO  lsp = 1, nvar
3516!
3517!--       Initialize
3518             vs = 0.0_wp
3519             vd_lu = 0.0_wp
3520             rs = 0.0_wp
3521             rb = 0.0_wp
3522             rc_tot = 0.0_wp
3523             IF ( spc_names(lsp) == 'PM10' )  THEN
3524                part_type = 1
3525!
3526!--             Sedimentation velocity:
3527                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3528                     particle_pars(ind_p_size, part_type), &
3529                     particle_pars(ind_p_slip, part_type), &
3530                     visc)
3531
3532                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3533                     vs, &
3534                     particle_pars(ind_p_size, part_type), &
3535                     particle_pars(ind_p_slip, part_type), &
3536                     nwet, temp_tmp, dens, visc, &
3537                     lup_dep,  &
3538                     r_aero_surf, ustar_surf )
3539
3540                bud_lup(lsp) = - conc_ijk(lsp) * &
3541                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3542
3543
3544             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3545                part_type = 2
3546!
3547!--             Sedimentation velocity:
3548                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3549                     particle_pars(ind_p_size, part_type), &
3550                     particle_pars(ind_p_slip, part_type), &
3551                     visc)
3552
3553                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3554                     vs, &
3555                     particle_pars(ind_p_size, part_type), &
3556                     particle_pars(ind_p_slip, part_type), &
3557                     nwet, temp_tmp, dens, visc, &
3558                     lup_dep, &
3559                     r_aero_surf, ustar_surf )
3560
3561                bud_lup(lsp) = - conc_ijk(lsp) * &
3562                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3563
3564             ELSE  !<GASES
3565!
3566!--             Read spc_name of current species for gas parameter
3567
3568                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
3569                   i_pspec = 0
3570                   DO  pspec = 1, nposp
3571                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3572                         i_pspec = pspec
3573                      END IF
3574                   ENDDO
3575
3576                ELSE
3577!
3578!--                For now species not deposited
3579                   CYCLE
3580                ENDIF
3581!
3582!--             Factor used for conversion from ppb to ug/m3 :
3583!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3584!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3585!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
3586!--             thus:
3587!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3588!--             Use density at lowest layer:
3589
3590                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3591!
3592!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3593                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3594                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
3595!
3596!--             Diffusivity for DEPAC relevant gases
3597!--             Use default value
3598                diffusivity            = 0.11e-4
3599!
3600!--             overwrite with known coefficients of diffusivity from Massman (1998)
3601                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3602                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3603                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3604                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3605                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3606                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3607                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3608!
3609!--             Get quasi-laminar boundary layer resistance rb:
3610                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland),   &
3611                     z0h_surf, ustar_surf, diffusivity, rb )
3612!
3613!--             Get rc_tot
3614                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,      &
3615                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep, 2,    &
3616                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,              &
3617                                         r_aero_surf , rb )
3618!
3619!--             Calculate budget
3620                IF ( rc_tot <= 0.0 )  THEN
3621                   bud_lup(lsp) = 0.0_wp
3622                ELSE
3623                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3624                   bud_lup(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3625                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3626                ENDIF
3627
3628             ENDIF
3629          ENDDO
3630       ENDIF
3631!
3632!--    Water
3633       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
3634!
3635!--       No vegetation on water:
3636          lai = 0.0_wp
3637          sai = 0.0_wp
3638          slinnfac = 1.0_wp
3639!
3640!--       Get vd
3641          DO  lsp = 1, nvar
3642!
3643!--          Initialize
3644             vs = 0.0_wp
3645             vd_lu = 0.0_wp
3646             rs = 0.0_wp
3647             rb = 0.0_wp
3648             rc_tot = 0.0_wp 
3649             IF ( spc_names(lsp) == 'PM10' )  THEN
3650                part_type = 1
3651!
3652!--             Sedimentation velocity:
3653                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3654                     particle_pars(ind_p_size, part_type), &
3655                     particle_pars(ind_p_slip, part_type), &
3656                     visc)
3657
3658                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3659                     vs, &
3660                     particle_pars(ind_p_size, part_type), &
3661                     particle_pars(ind_p_slip, part_type), &
3662                     nwet, temp_tmp, dens, visc, &
3663                     luw_dep, &
3664                     r_aero_surf, ustar_surf )
3665
3666                bud_luw(lsp) = - conc_ijk(lsp) * &
3667                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3668
3669             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3670                part_type = 2
3671!
3672!--             Sedimentation velocity:
3673                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3674                     particle_pars(ind_p_size, part_type), &
3675                     particle_pars(ind_p_slip, part_type), &
3676                     visc)
3677
3678                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3679                     vs, &
3680                     particle_pars(ind_p_size, part_type), &
3681                     particle_pars(ind_p_slip, part_type), &
3682                     nwet, temp_tmp, dens, visc, &
3683                     luw_dep, &
3684                     r_aero_surf, ustar_surf )
3685
3686                bud_luw(lsp) = - conc_ijk(lsp) * &
3687                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3688
3689             ELSE  !<GASES
3690!
3691!--             Read spc_name of current species for gas PARAMETER
3692
3693                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
3694                   i_pspec = 0
3695                   DO  pspec = 1, nposp
3696                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3697                         i_pspec = pspec
3698                      END IF
3699                   ENDDO
3700
3701                ELSE
3702!
3703!--                For now species not deposited
3704                   CYCLE
3705                ENDIF
3706!
3707!--             Factor used for conversion from ppb to ug/m3 :
3708!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3709!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3710!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
3711!--             thus:
3712!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3713!--             Use density at lowest layer:
3714
3715                ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3716!
3717!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3718!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3719                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
3720!
3721!--             Diffusivity for DEPAC relevant gases
3722!--             Use default value
3723                diffusivity            = 0.11e-4
3724!
3725!--             overwrite with known coefficients of diffusivity from Massman (1998)
3726                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3727                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3728                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3729                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3730                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3731                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3732                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3733!
3734!--             Get quasi-laminar boundary layer resistance rb:
3735                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland),  &
3736                     z0h_surf, ustar_surf, diffusivity, rb )
3737
3738!--             Get rc_tot
3739                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         & 
3740                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep, 2,    &
3741                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
3742                                         r_aero_surf , rb )
3743!
3744!--             Calculate budget
3745                IF ( rc_tot <= 0.0 )  THEN
3746
3747                   bud_luw(lsp) = 0.0_wp
3748
3749                ELSE
3750
3751                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3752
3753                   bud_luw(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3754                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3755                ENDIF
3756
3757             ENDIF
3758          ENDDO
3759       ENDIF
3760
3761
3762       bud = 0.0_wp
3763!
3764!--    Calculate overall budget for surface m and adapt concentration
3765       DO  lsp = 1, nspec
3766
3767          bud(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_luv(lsp) + &
3768               surf_lsm_h%frac(ind_pav_green,m) * bud_lup(lsp) + &
3769               surf_lsm_h%frac(ind_wat_win,m) * bud_luw(lsp)
3770!
3771!--       Compute new concentration:
3772          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
3773
3774          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, conc_ijk(lsp))
3775
3776       ENDDO
3777
3778    ENDIF
3779!
3780!-- For USM surfaces   
3781
3782    IF ( match_usm )  THEN
3783!
3784!--    Get surface element information at i,j:
3785       m = surf_usm_h%start_index(j,i)
3786       k = surf_usm_h%k(m)
3787!
3788!--    Get needed variables for surface element m
3789       ustar_surf  = surf_usm_h%us(m)
3790       z0h_surf    = surf_usm_h%z0h(m)
3791       r_aero_surf = surf_usm_h%r_a(m)
3792       solar_rad   = surf_usm_h%rad_sw_dir(m) + surf_usm_h%rad_sw_dif(m)
3793       lai = surf_usm_h%lai(m)
3794       sai = lai + 1
3795!
3796!--    For small grid spacing neglect R_a
3797       IF ( dzw(k) <= 1.0 )  THEN
3798          r_aero_surf = 0.0_wp
3799       ENDIF
3800!
3801!--    Initialize lu's
3802       luu_palm = 0
3803       luu_dep = 0
3804       lug_palm = 0
3805       lug_dep = 0
3806       lud_palm = 0
3807       lud_dep = 0
3808!
3809!--    Initialize budgets
3810       bud_luu  = 0.0_wp
3811       bud_lug = 0.0_wp
3812       bud_lud = 0.0_wp
3813!
3814!--    Get land use for i,j and assign to DEPAC lu
3815       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
3816!
3817!--       For green urban surfaces (e.g. green roofs
3818!--       assume LU short grass
3819          lug_palm = ind_luv_s_grass
3820          IF ( lug_palm == ind_luv_user )  THEN
3821             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
3822             CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
3823          ELSEIF ( lug_palm == ind_luv_b_soil )  THEN
3824             lug_dep = 9
3825          ELSEIF ( lug_palm == ind_luv_mixed_crops )  THEN
3826             lug_dep = 2
3827          ELSEIF ( lug_palm == ind_luv_s_grass )  THEN
3828             lug_dep = 1
3829          ELSEIF ( lug_palm == ind_luv_ev_needle_trees )  THEN
3830             lug_dep = 4
3831          ELSEIF ( lug_palm == ind_luv_de_needle_trees )  THEN
3832             lug_dep = 4
3833          ELSEIF ( lug_palm == ind_luv_ev_broad_trees )  THEN
3834             lug_dep = 12
3835          ELSEIF ( lug_palm == ind_luv_de_broad_trees )  THEN
3836             lug_dep = 5
3837          ELSEIF ( lug_palm == ind_luv_t_grass )  THEN
3838             lug_dep = 1
3839          ELSEIF ( lug_palm == ind_luv_desert )  THEN
3840             lug_dep = 9
3841          ELSEIF ( lug_palm == ind_luv_tundra  )  THEN
3842             lug_dep = 8
3843          ELSEIF ( lug_palm == ind_luv_irr_crops )  THEN
3844             lug_dep = 2
3845          ELSEIF ( lug_palm == ind_luv_semidesert )  THEN
3846             lug_dep = 8
3847          ELSEIF ( lug_palm == ind_luv_ice )  THEN
3848             lug_dep = 10
3849          ELSEIF ( lug_palm == ind_luv_marsh )  THEN
3850             lug_dep = 8
3851          ELSEIF ( lug_palm == ind_luv_ev_shrubs )  THEN
3852             lug_dep = 14
3853          ELSEIF ( lug_palm == ind_luv_de_shrubs  )  THEN
3854             lug_dep = 14
3855          ELSEIF ( lug_palm == ind_luv_mixed_forest )  THEN
3856             lug_dep = 4
3857          ELSEIF ( lug_palm == ind_luv_intrup_forest )  THEN
3858             lug_dep = 8     
3859          ENDIF
3860       ENDIF
3861
3862       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
3863!
3864!--       For walls in USM assume concrete walls/roofs,
3865!--       assumed LU class desert as also assumed for
3866!--       pavements in LSM
3867          luu_palm = ind_lup_conc
3868          IF ( luu_palm == ind_lup_user )  THEN
3869             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3870             CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 )
3871          ELSEIF ( luu_palm == ind_lup_asph_conc )  THEN
3872             luu_dep = 9
3873          ELSEIF ( luu_palm == ind_lup_asph )  THEN
3874             luu_dep = 9
3875          ELSEIF ( luu_palm ==  ind_lup_conc )  THEN
3876             luu_dep = 9
3877          ELSEIF ( luu_palm ==  ind_lup_sett )  THEN
3878             luu_dep = 9
3879          ELSEIF ( luu_palm == ind_lup_pav_stones )  THEN
3880             luu_dep = 9
3881          ELSEIF ( luu_palm == ind_lup_cobblest )  THEN
3882             luu_dep = 9       
3883          ELSEIF ( luu_palm == ind_lup_metal )  THEN
3884             luu_dep = 9
3885          ELSEIF ( luu_palm == ind_lup_wood )  THEN
3886             luu_dep = 9   
3887          ELSEIF ( luu_palm == ind_lup_gravel )  THEN
3888             luu_dep = 9
3889          ELSEIF ( luu_palm == ind_lup_f_gravel )  THEN
3890             luu_dep = 9
3891          ELSEIF ( luu_palm == ind_lup_pebblest )  THEN
3892             luu_dep = 9
3893          ELSEIF ( luu_palm == ind_lup_woodchips )  THEN
3894             luu_dep = 9
3895          ELSEIF ( luu_palm == ind_lup_tartan )  THEN
3896             luu_dep = 9
3897          ELSEIF ( luu_palm == ind_lup_art_turf )  THEN
3898             luu_dep = 9
3899          ELSEIF ( luu_palm == ind_lup_clay )  THEN
3900             luu_dep = 9
3901          ENDIF
3902       ENDIF
3903
3904       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
3905!
3906!--       For windows in USM assume metal as this is
3907!--       as close as we get, assumed LU class desert
3908!--       as also assumed for pavements in LSM
3909          lud_palm = ind_lup_metal     
3910          IF ( lud_palm == ind_lup_user )  THEN
3911             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3912             CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 )
3913          ELSEIF ( lud_palm == ind_lup_asph_conc )  THEN
3914             lud_dep = 9
3915          ELSEIF ( lud_palm == ind_lup_asph )  THEN
3916             lud_dep = 9
3917          ELSEIF ( lud_palm ==  ind_lup_conc )  THEN
3918             lud_dep = 9
3919          ELSEIF ( lud_palm ==  ind_lup_sett )  THEN
3920             lud_dep = 9
3921          ELSEIF ( lud_palm == ind_lup_pav_stones )  THEN
3922             lud_dep = 9
3923          ELSEIF ( lud_palm == ind_lup_cobblest )  THEN
3924             lud_dep = 9       
3925          ELSEIF ( lud_palm == ind_lup_metal )  THEN
3926             lud_dep = 9
3927          ELSEIF ( lud_palm == ind_lup_wood )  THEN
3928             lud_dep = 9   
3929          ELSEIF ( lud_palm == ind_lup_gravel )  THEN
3930             lud_dep = 9
3931          ELSEIF ( lud_palm == ind_lup_f_gravel )  THEN
3932             lud_dep = 9
3933          ELSEIF ( lud_palm == ind_lup_pebblest )  THEN
3934             lud_dep = 9
3935          ELSEIF ( lud_palm == ind_lup_woodchips )  THEN
3936             lud_dep = 9
3937          ELSEIF ( lud_palm == ind_lup_tartan )  THEN
3938             lud_dep = 9
3939          ELSEIF ( lud_palm == ind_lup_art_turf )  THEN
3940             lud_dep = 9
3941          ELSEIF ( lud_palm == ind_lup_clay )  THEN
3942             lud_dep = 9
3943          ENDIF
3944       ENDIF
3945!
3946!--    @TODO: Activate these lines as soon as new ebsolver branch is merged:
3947!--    Set wetness indicator to dry or wet for usm vegetation or pavement
3948       !IF ( surf_usm_h%c_liq(m) > 0 )  THEN
3949       !   nwet = 1
3950       !ELSE
3951       nwet = 0
3952       !ENDIF
3953!
3954!--    Compute length of time step
3955       IF ( call_chem_at_all_substeps )  THEN
3956          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
3957       ELSE
3958          dt_chem = dt_3d
3959       ENDIF
3960
3961       dh = dzw(k)
3962       inv_dh = 1.0_wp / dh
3963       dt_dh = dt_chem/dh
3964!
3965!--    Concentration at i,j,k
3966       DO  lsp = 1, nspec
3967          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
3968       ENDDO
3969!
3970!--    Temperature at i,j,k
3971       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
3972
3973       ts       = temp_tmp - 273.15  !< in degrees celcius
3974!
3975!--    Viscosity of air
3976       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
3977!
3978!--    Air density at k
3979       dens = rho_air_zw(k)
3980!
3981!--    Calculate relative humidity from specific humidity for DEPAC
3982       qv_tmp = MAX(q(k,j,i),0.0_wp)
3983       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
3984!
3985!--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
3986!--    for each surface fraction. Then derive overall budget taking into account the surface fractions.
3987!
3988!--    Walls/roofs
3989       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
3990!
3991!--       No vegetation on non-green walls:
3992          lai = 0.0_wp
3993          sai = 0.0_wp
3994
3995          slinnfac = 1.0_wp
3996!
3997!--       Get vd
3998          DO  lsp = 1, nvar
3999!
4000!--          Initialize
4001             vs = 0.0_wp
4002             vd_lu = 0.0_wp
4003             rs = 0.0_wp
4004             rb = 0.0_wp
4005             rc_tot = 0.0_wp
4006             IF (spc_names(lsp) == 'PM10' )  THEN
4007                part_type = 1
4008!
4009!--             Sedimentation velocity
4010                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4011                     particle_pars(ind_p_size, part_type), &
4012                     particle_pars(ind_p_slip, part_type), &
4013                     visc)
4014
4015                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4016                     vs, &
4017                     particle_pars(ind_p_size, part_type), &
4018                     particle_pars(ind_p_slip, part_type), &
4019                     nwet, temp_tmp, dens, visc, &
4020                     luu_dep,  &
4021                     r_aero_surf, ustar_surf )
4022
4023                bud_luu(lsp) = - conc_ijk(lsp) * &
4024                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4025
4026             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4027                part_type = 2
4028!
4029!--             Sedimentation velocity
4030                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4031                     particle_pars(ind_p_size, part_type), &
4032                     particle_pars(ind_p_slip, part_type), &
4033                     visc)
4034
4035                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4036                     vs, &
4037                     particle_pars(ind_p_size, part_type), &
4038                     particle_pars(ind_p_slip, part_type), &
4039                     nwet, temp_tmp, dens, visc, &
4040                     luu_dep , &
4041                     r_aero_surf, ustar_surf )
4042
4043                bud_luu(lsp) = - conc_ijk(lsp) * &
4044                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4045
4046             ELSE  !< GASES
4047!
4048!--             Read spc_name of current species for gas parameter
4049
4050                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4051                   i_pspec = 0
4052                   DO  pspec = 1, nposp
4053                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4054                         i_pspec = pspec
4055                      END IF
4056                   ENDDO
4057                ELSE
4058!
4059!--                For now species not deposited
4060                   CYCLE
4061                ENDIF
4062!
4063!--             Factor used for conversion from ppb to ug/m3 :
4064!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4065!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4066!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
4067!--             thus:
4068!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4069!--             Use density at k:
4070
4071                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
4072
4073                !
4074!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4075!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4076                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4077!
4078!--             Diffusivity for DEPAC relevant gases
4079!--             Use default value
4080                diffusivity            = 0.11e-4
4081!
4082!--             overwrite with known coefficients of diffusivity from Massman (1998)
4083                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4084                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4085                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4086                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4087                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4088                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4089                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4090!
4091!--             Get quasi-laminar boundary layer resistance rb:
4092                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland),   &
4093                     z0h_surf, ustar_surf, diffusivity, &
4094                     rb )
4095!
4096!--             Get rc_tot
4097                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,          &
4098                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, 2,     &
4099                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,   &
4100                                         r_aero_surf, rb )
4101!
4102!--             Calculate budget
4103                IF ( rc_tot <= 0.0 )  THEN
4104
4105                   bud_luu(lsp) = 0.0_wp
4106
4107                ELSE
4108
4109                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4110
4111                   bud_luu(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4112                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4113                ENDIF
4114
4115             ENDIF
4116          ENDDO
4117       ENDIF
4118!
4119!--    Green usm surfaces
4120       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
4121
4122          slinnfac = 1.0_wp
4123!
4124!--       Get vd
4125          DO  lsp = 1, nvar
4126!
4127!--          Initialize
4128             vs = 0.0_wp
4129             vd_lu = 0.0_wp
4130             rs = 0.0_wp
4131             rb = 0.0_wp
4132             rc_tot = 0.0_wp
4133             IF ( spc_names(lsp) == 'PM10' )  THEN
4134                part_type = 1
4135!
4136!--             Sedimentation velocity
4137                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4138                     particle_pars(ind_p_size, part_type), &
4139                     particle_pars(ind_p_slip, part_type), &
4140                     visc)
4141
4142                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4143                     vs, &
4144                     particle_pars(ind_p_size, part_type), &
4145                     particle_pars(ind_p_slip, part_type), &
4146                     nwet, temp_tmp, dens, visc, &
4147                     lug_dep,  &
4148                     r_aero_surf, ustar_surf )
4149
4150                bud_lug(lsp) = - conc_ijk(lsp) * &
4151                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4152
4153             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4154                part_type = 2
4155!
4156!--             Sedimentation velocity
4157                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4158                     particle_pars(ind_p_size, part_type), &
4159                     particle_pars(ind_p_slip, part_type), &
4160                     visc)
4161
4162                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4163                     vs, &
4164                     particle_pars(ind_p_size, part_type), &
4165                     particle_pars(ind_p_slip, part_type), &
4166                     nwet, temp_tmp, dens, visc, &
4167                     lug_dep, &
4168                     r_aero_surf, ustar_surf )
4169
4170                bud_lug(lsp) = - conc_ijk(lsp) * &
4171                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4172
4173             ELSE  !< GASES
4174!
4175!--             Read spc_name of current species for gas parameter
4176
4177                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4178                   i_pspec = 0
4179                   DO  pspec = 1, nposp
4180                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4181                         i_pspec = pspec
4182                      END IF
4183                   ENDDO
4184                ELSE
4185!
4186!--                For now species not deposited
4187                   CYCLE
4188                ENDIF
4189!
4190!--             Factor used for conversion from ppb to ug/m3 :
4191!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4192!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4193!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
4194!--             thus:
4195!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4196!--             Use density at k:
4197
4198                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
4199!
4200!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4201                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4202                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4203!
4204!--             Diffusivity for DEPAC relevant gases
4205!--             Use default value
4206                diffusivity            = 0.11e-4
4207!
4208!--             overwrite with known coefficients of diffusivity from Massman (1998)
4209                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4210                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4211                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4212                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4213                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4214                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4215                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4216!
4217!--             Get quasi-laminar boundary layer resistance rb:
4218                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland),    &
4219                     z0h_surf, ustar_surf, diffusivity, &
4220                     rb )
4221!
4222!--             Get rc_tot
4223                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,           &
4224                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, 2,      &
4225                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,    &
4226                                         r_aero_surf , rb )
4227!
4228!--             Calculate budget
4229                IF ( rc_tot <= 0.0 )  THEN
4230
4231                   bud_lug(lsp) = 0.0_wp
4232
4233                ELSE
4234
4235                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4236
4237                   bud_lug(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4238                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4239                ENDIF
4240
4241             ENDIF
4242          ENDDO
4243       ENDIF
4244!
4245!--    Windows
4246       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
4247!
4248!--       No vegetation on windows:
4249          lai = 0.0_wp
4250          sai = 0.0_wp
4251
4252          slinnfac = 1.0_wp
4253!
4254!--       Get vd
4255          DO  lsp = 1, nvar
4256!
4257!--          Initialize
4258             vs = 0.0_wp
4259             vd_lu = 0.0_wp
4260             rs = 0.0_wp
4261             rb = 0.0_wp
4262             rc_tot = 0.0_wp 
4263             IF ( spc_names(lsp) == 'PM10' )  THEN
4264                part_type = 1
4265!
4266!--             Sedimentation velocity
4267                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4268                     particle_pars(ind_p_size, part_type), &
4269                     particle_pars(ind_p_slip, part_type), &
4270                     visc)
4271
4272                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
4273                     particle_pars(ind_p_size, part_type), &
4274                     particle_pars(ind_p_slip, part_type), &
4275                     nwet, temp_tmp, dens, visc,              &
4276                     lud_dep, r_aero_surf, ustar_surf )
4277
4278                bud_lud(lsp) = - conc_ijk(lsp) * &
4279                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4280
4281             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4282                part_type = 2
4283!
4284!--             Sedimentation velocity
4285                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4286                     particle_pars(ind_p_size, part_type), &
4287                     particle_pars(ind_p_slip, part_type), &
4288                     visc)
4289
4290                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
4291                     particle_pars(ind_p_size, part_type), &
4292                     particle_pars(ind_p_slip, part_type), &
4293                     nwet, temp_tmp, dens, visc, &
4294                     lud_dep, &
4295                     r_aero_surf, ustar_surf )
4296
4297                bud_lud(lsp) = - conc_ijk(lsp) * &
4298                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4299
4300             ELSE  !< GASES
4301!
4302!--             Read spc_name of current species for gas PARAMETER
4303
4304                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4305                   i_pspec = 0
4306                   DO  pspec = 1, nposp
4307                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4308                         i_pspec = pspec
4309                      END IF
4310                   ENDDO
4311                ELSE
4312!
4313!--                For now species not deposited
4314                   CYCLE
4315                ENDIF
4316!
4317!--             Factor used for conversion from ppb to ug/m3 :
4318!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4319!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4320!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
4321!--             thus:
4322!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4323!--             Use density at k:
4324
4325                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
4326!
4327!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4328!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4329                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4330!
4331!--             Diffusivity for DEPAC relevant gases
4332!--             Use default value
4333                diffusivity = 0.11e-4
4334!
4335!--             overwrite with known coefficients of diffusivity from Massman (1998)
4336                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4337                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4338                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4339                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4340                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4341                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4342                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4343!
4344!--             Get quasi-laminar boundary layer resistance rb:
4345                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland),   &
4346                     z0h_surf, ustar_surf, diffusivity, rb )
4347!
4348!--             Get rc_tot
4349                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
4350                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, 2,    &
4351                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
4352                                         r_aero_surf , rb )
4353!
4354!--             Calculate budget
4355                IF ( rc_tot <= 0.0 )  THEN
4356
4357                   bud_lud(lsp) = 0.0_wp
4358
4359                ELSE
4360
4361                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4362
4363                   bud_lud(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4364                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4365                ENDIF
4366
4367             ENDIF
4368          ENDDO
4369       ENDIF
4370
4371
4372       bud = 0.0_wp
4373!
4374!--    Calculate overall budget for surface m and adapt concentration
4375       DO  lsp = 1, nspec
4376
4377
4378          bud(lsp) = surf_usm_h%frac(ind_veg_wall,m) * bud_luu(lsp) + &
4379               surf_usm_h%frac(ind_pav_green,m) * bud_lug(lsp) + &
4380               surf_usm_h%frac(ind_wat_win,m) * bud_lud(lsp)
4381!
4382!--       Compute new concentration
4383          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
4384
4385          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
4386
4387       ENDDO
4388
4389    ENDIF
4390
4391
4392 END SUBROUTINE chem_depo
4393
4394
4395!------------------------------------------------------------------------------!
4396! Description:
4397! ------------
4398!> Subroutine to compute total canopy (or surface) resistance Rc for gases
4399!>
4400!> DEPAC:
4401!> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
4402!> module of the LOTOS-EUROS model (Manders et al., 2017)
4403!>
4404!> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
4405!> van Zanten et al., 2010.
4406!------------------------------------------------------------------------------!
4407 SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi,    &
4408      rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,  &
4409      ra, rb ) 
4410!
4411!--   Some of depac arguments are OPTIONAL:
4412!--    A. compute Rc_tot without compensation points (ccomp_tot will be zero):
4413!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
4414!--    B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
4415!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4416!--                c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
4417!--
4418!--    C. compute effective Rc based on compensation points (used for OPS):
4419!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4420!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4421!--               ra, rb, rc_eff)
4422!--    X1. Extra (OPTIONAL) output variables:
4423!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4424!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4425!--               ra, rb, rc_eff, &
4426!--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
4427!--    X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
4428!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4429!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4430!--               ra, rb, rc_eff, &
4431!--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
4432!--               calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
4433
4434
4435    CHARACTER(LEN=*), INTENT(IN) ::  compnam         !< component name
4436                                                     !< 'HNO3','NO','NO2','O3','SO2','NH3'
4437    INTEGER(iwp), INTENT(IN) ::  day_of_year         !< day of year, 1 ... 365 (366)
4438    INTEGER(iwp), INTENT(IN) ::  nwet                !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4439    INTEGER(iwp), INTENT(IN) ::  lu                  !< land use type, lu = 1,...,nlu
4440    INTEGER(iwp), INTENT(IN) ::  iratns              !< index for NH3/SO2 ratio used for SO2:
4441                                                     !< iratns = 1: low NH3/SO2
4442                                                     !< iratns = 2: high NH3/SO2
4443                                                     !< iratns = 3: very low NH3/SO2
4444    REAL(wp), INTENT(IN) ::  lat                     !< latitude Northern hemisphere (degrees) (S. hemisphere not possible)
4445    REAL(wp), INTENT(IN) ::  t                       !< temperature (C)
4446    REAL(wp), INTENT(IN) ::  ust                     !< friction velocity (m/s)
4447    REAL(wp), INTENT(IN) ::  solar_rad               !< solar radiation, dirict+diffuse (W/m2)
4448    REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
4449    REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
4450    REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
4451    REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and stems)
4452    REAL(wp), INTENT(IN) ::  diffusivity             !< diffusivity
4453    REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
4454    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3           !< actual atmospheric concentration (ug/m3), in DEPAC=Catm
4455    REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
4456    REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
4457
4458    REAL(wp), INTENT(OUT) ::  rc_tot                 !< total canopy resistance Rc (s/m)
4459    REAL(wp), INTENT(OUT) ::  ccomp_tot              !< total compensation point (ug/m3)
4460!                                                     !< [= 0 for species that don't have a compensation
4461!-- Local variables:
4462!
4463!-- Component number taken from component name, paramteres matched with include files
4464    INTEGER(iwp) ::  icmp
4465!
4466!-- Component numbers:
4467    INTEGER(iwp), PARAMETER ::  icmp_o3   = 1
4468    INTEGER(iwp), PARAMETER ::  icmp_so2  = 2
4469    INTEGER(iwp), PARAMETER ::  icmp_no2  = 3
4470    INTEGER(iwp), PARAMETER ::  icmp_no   = 4
4471    INTEGER(iwp), PARAMETER ::  icmp_nh3  = 5
4472    INTEGER(iwp), PARAMETER ::  icmp_co   = 6
4473    INTEGER(iwp), PARAMETER ::  icmp_no3  = 7
4474    INTEGER(iwp), PARAMETER ::  icmp_hno3 = 8
4475    INTEGER(iwp), PARAMETER ::  icmp_n2o5 = 9
4476    INTEGER(iwp), PARAMETER ::  icmp_h2o2 = 10
4477
4478    LOGICAL ::  ready                                !< Rc has been set:
4479    !< = 1 -> constant Rc
4480    !< = 2 -> temperature dependent Rc
4481!
4482!-- Vegetation indicators:
4483    LOGICAL ::  lai_present                          !< leaves are present for current land use type
4484    LOGICAL ::  sai_present                          !< vegetation is present for current land use type
4485
4486!    REAL(wp) ::  laimax                              !< maximum leaf area index (-)
4487    REAL(wp) ::  gw                                  !< external leaf conductance (m/s)
4488    REAL(wp) ::  gstom                               !< stomatal conductance (m/s)
4489    REAL(wp) ::  gsoil_eff                           !< effective soil conductance (m/s)
4490    REAL(wp) ::  gc_tot                              !< total canopy conductance (m/s)
4491    REAL(wp) ::  cw                                  !< external leaf surface compensation point (ug/m3)
4492    REAL(wp) ::  cstom                               !< stomatal compensation point (ug/m3)
4493    REAL(wp) ::  csoil                               !< soil compensation point (ug/m3)
4494!
4495!-- Next statement is just to avoid compiler warning about unused variable
4496    IF ( day_of_year == 0  .OR.  ( conc_ijk_ugm3 + lat + ra + rb ) > 0.0_wp )  CONTINUE
4497!
4498!-- Define component number
4499    SELECT CASE ( TRIM( compnam ) )
4500
4501    CASE ( 'O3', 'o3' )
4502       icmp = icmp_o3
4503
4504    CASE ( 'SO2', 'so2' )
4505       icmp = icmp_so2
4506
4507    CASE ( 'NO2', 'no2' )
4508       icmp = icmp_no2
4509
4510    CASE ( 'NO', 'no' )
4511       icmp = icmp_no
4512
4513    CASE ( 'NH3', 'nh3' )
4514       icmp = icmp_nh3
4515
4516    CASE ( 'CO', 'co' )
4517       icmp = icmp_co
4518
4519    CASE ( 'NO3', 'no3' )
4520       icmp = icmp_no3
4521
4522    CASE ( 'HNO3', 'hno3' )
4523       icmp = icmp_hno3
4524
4525    CASE ( 'N2O5', 'n2o5' )
4526       icmp = icmp_n2o5
4527
4528    CASE ( 'H2O2', 'h2o2' )
4529       icmp = icmp_h2o2
4530
4531    CASE default
4532!
4533!--    Component not part of DEPAC --> not deposited
4534       RETURN
4535
4536    END SELECT
4537
4538!
4539!-- Inititalize
4540    gw        = 0.0_wp
4541    gstom     = 0.0_wp
4542    gsoil_eff = 0.0_wp
4543    gc_tot    = 0.0_wp
4544    cw        = 0.0_wp
4545    cstom     = 0.0_wp
4546    csoil     = 0.0_wp
4547!
4548!-- Check whether vegetation is present:
4549    lai_present = ( lai > 0.0 )
4550    sai_present = ( sai > 0.0 )
4551!
4552!-- Set Rc (i.e. rc_tot) in special cases:
4553    CALL rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
4554!
4555!-- If Rc is not set:
4556    IF ( .NOT. ready ) then
4557!
4558!--    External conductance:
4559       CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )         
4560!
4561!--    Stomatal conductance:
4562       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p )
4563!
4564!--    Effective soil conductance:
4565       CALL rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff )
4566!
4567!--    Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
4568       CALL rc_rctot( gstom, gsoil_eff, gw, gc_tot, rc_tot )
4569!
4570!--    Needed to include compensation point for NH3
4571!--    Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
4572!--    CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
4573!--          lai_present, sai_present, &
4574!--          ccomp_tot, &
4575!--          conc_ijk_ugm3=conc_ijk_ugm3,c_ave_prev_nh3=c_ave_prev_nh3, &
4576!--          c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
4577!--          tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
4578!
4579!--    Effective Rc based on compensation points:
4580!--        IF ( present(rc_eff) ) then
4581!--          check on required arguments:
4582!--           IF ( (.not. present(conc_ijk_ugm3)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
4583!--              stop 'output argument rc_eff requires input arguments conc_ijk_ugm3, ra and rb'
4584!--           END IF
4585!
4586!--       Compute rc_eff :
4587       !      CALL rc_comp_point_rc_eff(ccomp_tot,conc_ijk_ugm3,ra,rb,rc_tot,rc_eff)
4588       !   ENDIF
4589    ENDIF
4590
4591 END SUBROUTINE drydepos_gas_depac
4592
4593
4594!------------------------------------------------------------------------------!
4595! Description:
4596! ------------
4597!> Subroutine to compute total canopy resistance in special cases
4598!------------------------------------------------------------------------------!
4599 SUBROUTINE rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
4600
4601   
4602    CHARACTER(LEN=*), INTENT(IN)  ::  compnam     !< component name
4603
4604    INTEGER(iwp), INTENT(IN)  ::  icmp            !< component index     
4605    INTEGER(iwp), INTENT(IN)  ::  lu              !< land use type, lu = 1,...,nlu
4606    INTEGER(iwp), INTENT(IN)  ::  nwet            !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4607
4608    REAL(wp), INTENT(IN)  ::  t                   !< temperature (C)
4609
4610    REAL(wp), INTENT(OUT) ::  rc_tot             !< total canopy resistance Rc (s/m)
4611    REAL(wp), INTENT(OUT) ::  ccomp_tot          !< total compensation point (ug/m3)
4612
4613    LOGICAL, INTENT(OUT) ::  ready               !< Rc has been set
4614                                                 !< = 1 -> constant Rc
4615!
4616!-- Next line is to avoid compiler warning about unused variable
4617    IF ( icmp == 0 )  CONTINUE
4618!
4619!-- rc_tot is not yet set:
4620    ready = .FALSE.
4621!
4622!-- Default compensation point in special CASEs = 0:
4623    ccomp_tot = 0.0_wp
4624
4625    SELECT CASE( TRIM( compnam ) )
4626    CASE( 'HNO3', 'N2O5', 'NO3', 'H2O2' )
4627!
4628!--    No separate resistances for HNO3; just one total canopy resistance:
4629       IF ( t < -5.0_wp .AND. nwet == 9 )  THEN
4630!
4631!--       T < 5 C and snow:
4632          rc_tot = 50.0_wp
4633       ELSE
4634!
4635!--       all other circumstances:
4636          rc_tot = 10.0_wp
4637       ENDIF
4638       ready = .TRUE.
4639
4640    CASE( 'NO', 'CO' )
4641       IF ( lu == ilu_water_sea .OR. lu == ilu_water_inland )  THEN       ! water
4642          rc_tot = 2000.0_wp
4643          ready = .TRUE.
4644       ELSEIF ( nwet == 1 )  THEN  !< wet
4645          rc_tot = 2000.0_wp
4646          ready = .TRUE.
4647       ENDIF
4648    CASE( 'NO2', 'O3', 'SO2', 'NH3' )
4649!
4650!--    snow surface:
4651       IF ( nwet == 9 )  THEN
4652!
4653!--       To be activated when snow is implemented
4654          !CALL rc_snow(ipar_snow(icmp),t,rc_tot)
4655          ready = .TRUE.
4656       ENDIF
4657    CASE default
4658       message_string = 'Component '// TRIM( compnam ) // ' not supported'
4659       CALL message( 'rc_special', 'CM0457', 1, 2, 0, 6, 0 )
4660    END SELECT
4661
4662 END SUBROUTINE rc_special
4663
4664
4665!------------------------------------------------------------------------------!
4666! Description:
4667! ------------
4668!> Subroutine to compute external conductance
4669!------------------------------------------------------------------------------!
4670 SUBROUTINE rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai, gw )
4671
4672!
4673!-- Input/output variables:
4674    CHARACTER(LEN=*), INTENT(IN) ::  compnam      !< component name
4675
4676    INTEGER(iwp), INTENT(IN) ::  nwet             !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4677    INTEGER(iwp), INTENT(IN) ::  iratns           !< index for NH3/SO2 ratio;
4678                                                  !< iratns = 1: low NH3/SO2
4679                                                  !< iratns = 2: high NH3/SO2
4680                                                  !< iratns = 3: very low NH3/SO2
4681    LOGICAL, INTENT(IN) ::  sai_present
4682
4683    REAL(wp), INTENT(IN) ::  t                    !< temperature (C)
4684    REAL(wp), INTENT(IN) ::  rh                   !< relative humidity (%)
4685    REAL(wp), INTENT(IN) ::  sai                  !< one-sided leaf area index (-)
4686
4687    REAL(wp), INTENT(OUT) ::  gw                  !< external leaf conductance (m/s)
4688
4689    SELECT CASE( TRIM( compnam ) )
4690
4691    CASE( 'NO2' )
4692       CALL rw_constant( 2000.0_wp, sai_present, gw )
4693
4694    CASE( 'NO', 'CO' )
4695       CALL rw_constant( -9999.0_wp, sai_present, gw )   !< see Erisman et al, 1994 section 3.2.3
4696
4697    CASE( 'O3' )
4698       CALL rw_constant( 2500.0_wp, sai_present, gw )
4699
4700    CASE( 'SO2' )
4701       CALL rw_so2( t, nwet, rh, iratns, sai_present, gw )
4702
4703    CASE( 'NH3' )
4704       CALL rw_nh3_sutton( t, rh, sai_present, gw )
4705!
4706!--    conversion from leaf resistance to canopy resistance by multiplying with sai:
4707       gw = sai * gw
4708
4709    CASE default
4710       message_string = 'Component '// TRIM( compnam ) // ' not supported'
4711       CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 )
4712    END SELECT
4713
4714 END SUBROUTINE rc_gw
4715
4716
4717!------------------------------------------------------------------------------!
4718! Description:
4719! ------------
4720!> Subroutine to compute external leaf conductance for SO2
4721!------------------------------------------------------------------------------!
4722 SUBROUTINE rw_so2( t, nwet, rh, iratns, sai_present, gw )
4723
4724!
4725!-- Input/output variables:
4726    INTEGER(iwp), INTENT(IN) ::  nwet        !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4727    INTEGER(iwp), INTENT(IN) ::  iratns      !< index for NH3/SO2 ratio:
4728                                             !< iratns = 1: low NH3/SO2
4729                                             !< iratns = 2: high NH3/SO2
4730                                             !< iratns = 3: very low NH3/SO2
4731    LOGICAL, INTENT(IN) ::  sai_present
4732
4733    REAL(wp), INTENT(IN) ::  t               !< temperature (C)
4734    REAL(wp), INTENT(IN) ::  rh              !< relative humidity (%)   
4735
4736    REAL(wp), INTENT(OUT) ::  gw             !< external leaf conductance (m/s)
4737!
4738!-- Local variables:
4739    REAL(wp) ::  rw                          !< external leaf resistance (s/m)
4740!
4741!-- Check if vegetation present:
4742    IF ( sai_present )  THEN
4743
4744       IF ( nwet == 0 )  THEN
4745!
4746!--   ------------------------
4747!--         dry surface
4748!--   ------------------------
4749!--         T > -1 C
4750          IF ( t > -1.0_wp )  THEN
4751             IF ( rh < 81.3_wp )  THEN
4752                rw = 25000.0_wp * exp( -0.0693_wp * rh )
4753             ELSE
4754                rw = 0.58e12 * exp( -0.278_wp * rh ) + 10.0_wp
4755             ENDIF
4756          ELSE
4757             ! -5 C < T <= -1 C
4758             IF ( t > -5.0_wp )  THEN
4759                rw = 200.0_wp
4760             ELSE
4761                ! T <= -5 C
4762                rw = 500.0_wp
4763             ENDIF
4764          ENDIF
4765       ELSE
4766!
4767!--   ------------------------
4768!--         wet surface
4769!--   ------------------------
4770          rw = 10.0_wp !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10
4771       ENDIF
4772!
4773!--    very low NH3/SO2 ratio:
4774       IF ( iratns == iratns_very_low ) rw = rw + 50.0_wp
4775!
4776!--      Conductance:
4777       gw = 1.0_wp / rw
4778    ELSE
4779!
4780!--      no vegetation:
4781       gw = 0.0_wp
4782    ENDIF
4783
4784 END SUBROUTINE rw_so2
4785
4786
4787!------------------------------------------------------------------------------!
4788! Description:
4789! ------------
4790!> Subroutine to compute external leaf conductance for NH3,
4791!>                  following Sutton & Fowler, 1993
4792!------------------------------------------------------------------------------!
4793 SUBROUTINE rw_nh3_sutton( tsurf, rh,sai_present, gw )
4794
4795!
4796!-- Input/output variables:
4797    LOGICAL, INTENT(IN) ::  sai_present
4798
4799    REAL(wp), INTENT(IN) ::  tsurf          !< surface temperature (C)
4800    REAL(wp), INTENT(IN) ::  rh             !< relative humidity (%)
4801
4802    REAL(wp), INTENT(OUT) ::  gw            !< external leaf conductance (m/s)
4803!
4804!-- Local variables:
4805    REAL(wp) ::  rw                         !< external leaf resistance (s/m)
4806    REAL(wp) ::  sai_grass_haarweg          !< surface area index at experimental site Haarweg
4807!
4808!-- Fix sai_grass at value valid for Haarweg data for which gamma_w parametrization is derived
4809    sai_grass_haarweg = 3.5_wp
4810!
4811!-- Calculation rw:
4812!--                    100 - rh
4813!--    rw = 2.0 * exp(----------)
4814!--                      12
4815
4816    IF ( sai_present )  THEN
4817!
4818!--    External resistance according to Sutton & Fowler, 1993
4819       rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp )
4820       rw = sai_grass_haarweg * rw
4821!
4822!--    Frozen soil (from Depac v1):
4823       IF ( tsurf < 0.0_wp ) rw = 200.0_wp
4824!
4825!--    Conductance:
4826       gw = 1.0_wp / rw
4827    ELSE
4828       ! no vegetation:
4829       gw = 0.0_wp
4830    ENDIF
4831
4832 END SUBROUTINE rw_nh3_sutton
4833
4834
4835!------------------------------------------------------------------------------!
4836! Description:
4837! ------------
4838!> Subroutine to compute external leaf conductance
4839!------------------------------------------------------------------------------!
4840 SUBROUTINE rw_constant( rw_val, sai_present, gw )
4841
4842!
4843!-- Input/output variables:
4844    LOGICAL, INTENT(IN) ::  sai_present
4845
4846    REAL(wp), INTENT(IN) ::  rw_val       !< constant value of Rw   
4847
4848    REAL(wp), INTENT(OUT) ::  gw          !< wernal leaf conductance (m/s)
4849!
4850!-- Compute conductance:
4851    IF ( sai_present .AND. .NOT.missing(rw_val) )  THEN
4852       gw = 1.0_wp / rw_val
4853    ELSE
4854       gw = 0.0_wp
4855    ENDIF
4856
4857 END SUBROUTINE rw_constant
4858
4859
4860!------------------------------------------------------------------------------!
4861! Description:
4862! ------------
4863!> Subroutine to compute stomatal conductance
4864!------------------------------------------------------------------------------!
4865 SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p ) 
4866
4867!
4868!-- input/output variables:
4869    CHARACTER(LEN=*), INTENT(IN) ::  compnam       !< component name
4870
4871    INTEGER(iwp), INTENT(IN) ::  icmp              !< component index
4872    INTEGER(iwp), INTENT(IN) ::  lu                !< land use type , lu = 1,...,nlu
4873
4874    LOGICAL, INTENT(IN) ::  lai_present
4875
4876    REAL(wp), INTENT(IN) ::  lai                   !< one-sided leaf area index
4877    REAL(wp), INTENT(IN) ::  solar_rad             !< solar radiation, dirict+diffuse (W/m2)
4878    REAL(wp), INTENT(IN) ::  sinphi                !< sin of solar elevation angle
4879    REAL(wp), INTENT(IN) ::  t                     !< temperature (C)
4880    REAL(wp), INTENT(IN) ::  rh                    !< relative humidity (%)
4881    REAL(wp), INTENT(IN) ::  diffusivity           !< diffusion coefficient of the gas involved
4882
4883    REAL(wp), OPTIONAL,INTENT(IN) :: p             !< pressure (Pa)
4884
4885    REAL(wp), INTENT(OUT) ::  gstom                !< stomatal conductance (m/s)
4886!
4887!-- Local variables
4888    REAL(wp) ::  vpd                               !< vapour pressure deficit (kPa)
4889
4890    REAL(wp), PARAMETER ::  dO3 = 0.13e-4          !< diffusion coefficient of ozon (m2/s)
4891!
4892!-- Next line is to avoid compiler warning about unused variables
4893    IF ( icmp == 0 )  CONTINUE
4894
4895    SELECT CASE( TRIM( compnam ) )
4896
4897    CASE( 'NO', 'CO' )
4898!
4899!--    For no stomatal uptake is neglected:
4900       gstom = 0.0_wp
4901
4902    CASE( 'NO2', 'O3', 'SO2', 'NH3' )
4903!
4904!--    if vegetation present:
4905       IF ( lai_present )  THEN
4906
4907          IF ( solar_rad > 0.0_wp )  THEN
4908             CALL rc_get_vpd( t, rh, vpd )
4909             CALL rc_gstom_emb( lu, solar_rad, t, vpd, lai_present, lai, sinphi, gstom, p )
4910             gstom = gstom * diffusivity / dO3       !< Gstom of Emberson is derived for ozone
4911          ELSE
4912             gstom = 0.0_wp
4913          ENDIF
4914       ELSE
4915!
4916!--       no vegetation; zero conductance (infinite resistance):
4917          gstom = 0.0_wp
4918       ENDIF
4919
4920    CASE default
4921       message_string = 'Component '// TRIM( compnam ) // ' not supported'
4922       CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 )
4923    END SELECT
4924
4925 END SUBROUTINE rc_gstom
4926
4927
4928!------------------------------------------------------------------------------!
4929! Description:
4930! ------------
4931!> Subroutine to compute stomatal conductance according to Emberson
4932!------------------------------------------------------------------------------!
4933 SUBROUTINE rc_gstom_emb( lu, solar_rad, T, vpd, lai_present, lai, sinp, Gsto, p )
4934!
4935!>  History
4936!>   Original code from Lotos-Euros, TNO, M. Schaap
4937!>   2009-08, M.C. van Zanten, Rivm
4938!>     Updated and extended.
4939!>   2009-09, Arjo Segers, TNO
4940!>     Limitted temperature influence to range to avoid
4941!>     floating point exceptions.
4942
4943!> Method
4944
4945!>   Code based on Emberson et al, 2000, Env. Poll., 403-413
4946!>   Notation conform Unified EMEP Model Description Part 1, ch 8
4947!
4948!>   In the calculation of f_light the modification of L. Zhang 2001, AE to the PARshade and PARsun
4949!>   parametrizations of Norman 1982 are applied
4950!>   f_phen and f_SWP are set to 1
4951!
4952!>   Land use types DEPAC versus Emberson (Table 5.1, EMEP model description)
4953!>   DEPAC                     Emberson
4954!>     1 = grass                 GR = grassland
4955!>     2 = arable land           TC = temperate crops ( lai according to RC = rootcrops)
4956!>     3 = permanent crops       TC = temperate crops ( lai according to RC = rootcrops)
4957!>     4 = coniferous forest     CF = tempareate/boREAL(wp) coniferous forest
4958!>     5 = deciduous forest      DF = temperate/boREAL(wp) deciduous forest
4959!>     6 = water                 W  = water
4960!>     7 = urban                 U  = urban
4961!>     8 = other                 GR = grassland
4962!>     9 = desert                DE = desert
4963!
4964!-- Emberson specific declarations
4965!
4966!-- Input/output variables:
4967    INTEGER(iwp), INTENT(IN) ::  lu             !< land use type, lu = 1,...,nlu
4968
4969    LOGICAL, INTENT(IN) ::  lai_present
4970
4971    REAL(wp), INTENT(IN) ::  solar_rad          !< solar radiation, dirict+diffuse (W/m2)
4972    REAL(wp), INTENT(IN) ::  t                  !< temperature (C)
4973    REAL(wp), INTENT(IN) ::  vpd                !< vapour pressure deficit (kPa)
4974
4975    REAL(wp), INTENT(IN) ::  lai                !< one-sided leaf area index
4976    REAL(wp), INTENT(IN) ::  sinp               !< sin of solar elevation angle
4977
4978    REAL(wp), OPTIONAL, INTENT(IN) ::  p        !< pressure (Pa)
4979
4980    REAL(wp), INTENT(OUT) ::  gsto              !< stomatal conductance (m/s)
4981!
4982!-- Local variables:
4983    REAL(wp) ::  f_light
4984    REAL(wp) ::  f_phen
4985    REAL(wp) ::  f_temp
4986    REAL(wp) ::  f_vpd
4987    REAL(wp) ::  f_swp
4988    REAL(wp) ::  bt
4989    REAL(wp) ::  f_env
4990    REAL(wp) ::  pardir
4991    REAL(wp) ::  pardiff
4992    REAL(wp) ::  parshade
4993    REAL(wp) ::  parsun
4994    REAL(wp) ::  laisun
4995    REAL(wp) ::  laishade
4996    REAL(wp) ::  sinphi
4997    REAL(wp) ::  pres
4998    REAL(wp), PARAMETER ::  p_sealevel = 1.01325e05    !< Pa
4999!
5000!-- Check whether vegetation is present:
5001    IF ( lai_present )  THEN
5002
5003       ! calculation of correction factors for stomatal conductance
5004       IF ( sinp <= 0.0_wp )  THEN 
5005          sinphi = 0.0001_wp
5006       ELSE
5007          sinphi = sinp
5008       END IF
5009!
5010!--    ratio between actual and sea-level pressure is used
5011!--    to correct for height in the computation of par;
5012!--    should not exceed sea-level pressure therefore ...
5013       IF (  present(p) )  THEN
5014          pres = min( p, p_sealevel )
5015       ELSE
5016          pres = p_sealevel
5017       ENDIF
5018!
5019!--    direct and diffuse par, Photoactive (=visible) radiation:
5020       CALL par_dir_diff( solar_rad, sinphi, pres, p_sealevel, pardir, pardiff )
5021!
5022!--    par for shaded leaves (canopy averaged):
5023       parshade = pardiff * exp( -0.5 * lai**0.7 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )     !< Norman,1982
5024       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
5025          parshade = pardiff * exp( -0.5 * lai**0.8 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )  !< Zhang et al., 2001
5026       END IF
5027!
5028!--    par for sunlit leaves (canopy averaged):
5029!--    alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
5030       parsun = pardir * 0.5/sinphi + parshade             !< Norman, 1982
5031       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
5032          parsun = pardir**0.8 * 0.5 / sinphi + parshade   !< Zhang et al., 2001
5033       END IF
5034!
5035!--    leaf area index for sunlit and shaded leaves:
5036       IF ( sinphi > 0 )  THEN
5037          laisun = 2 * sinphi * ( 1 - exp( -0.5 * lai / sinphi ) )
5038          laishade = lai - laisun
5039       ELSE
5040          laisun = 0
5041          laishade = lai
5042       END IF
5043
5044       f_light = ( laisun * ( 1 - exp( -1.0_wp * alpha(lu) * parsun ) ) + &
5045            laishade * ( 1 - exp( -1.0_wp * alpha(lu) * parshade ) ) ) / lai
5046
5047       f_light = MAX(f_light,f_min(lu))
5048!
5049!--    temperature influence; only non-zero within range [tmin,tmax]:
5050       IF ( ( tmin(lu) < t ) .AND. ( t < tmax(lu) ) )  THEN
5051          bt = ( tmax(lu) - topt(lu) ) / ( topt(lu) - tmin(lu) )
5052          f_temp = ( ( t - tmin(lu) ) / ( topt(lu) - tmin(lu) ) ) * ( ( tmax(lu) - t ) / ( tmax(lu) - topt(lu) ) )**bt
5053       ELSE
5054          f_temp = 0.0_wp
5055       END IF
5056       f_temp = MAX( f_temp, f_min(lu) )
5057!
5058!--      vapour pressure deficit influence
5059       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) ) )
5060       f_vpd = MAX( f_vpd, f_min(lu) )
5061
5062       f_swp = 1.0_wp
5063!
5064!--      influence of phenology on stom. conductance
5065!--      ignored for now in DEPAC since influence of f_phen on lu classes in use is negligible.
5066!--      When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore
5067       f_phen = 1.0_wp
5068!
5069!--      evaluate total stomatal conductance
5070       f_env = f_temp * f_vpd * f_swp
5071       f_env = MAX( f_env,f_min(lu) )
5072       gsto = g_max(lu) * f_light * f_phen * f_env
5073!
5074!--      gstom expressed per m2 leafarea;
5075!--      this is converted with lai to m2 surface.
5076       gsto = lai * gsto    ! in m/s
5077
5078    ELSE
5079       gsto = 0.0_wp
5080    ENDIF
5081
5082 END SUBROUTINE rc_gstom_emb
5083
5084
5085 !-------------------------------------------------------------------
5086 !> par_dir_diff
5087 !>     Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and
5088 !>     diffuse, visible and near-infrared components. Agric. Forest Meteorol.
5089 !>     34, 205-213.
5090 !>     From a SUBROUTINE obtained from Leiming Zhang,
5091 !>     Meteorological Service of Canada
5092 !>     Leiming uses solar irradiance. This should be equal to global radiation and
5093 !>     Willem Asman set it to global radiation (here defined as solar radiation, dirict+diffuse)
5094 !>
5095 !>     @todo Check/connect/replace with radiation_model_mod variables   
5096 !-------------------------------------------------------------------
5097 SUBROUTINE par_dir_diff( solar_rad, sinphi, pres, pres_0, par_dir, par_diff )
5098
5099
5100    REAL(wp), INTENT(IN) ::  solar_rad       !< solar radiation, dirict+diffuse (W m-2)
5101    REAL(wp), INTENT(IN) ::  sinphi          !< sine of the solar elevation
5102    REAL(wp), INTENT(IN) ::  pres            !< actual pressure (to correct for height) (Pa)
5103    REAL(wp), INTENT(IN) ::  pres_0          !< pressure at sea level (Pa)
5104
5105    REAL(wp), INTENT(OUT) ::  par_dir        !< par direct : visible (photoactive) direct beam radiation (W m-2)
5106    REAL(wp), INTENT(OUT) ::  par_diff       !< par diffuse: visible (photoactive) diffuse radiation (W m-2)
5107
5108
5109    REAL(wp) ::  sv                          !< total visible radiation
5110    REAL(wp) ::  fv                          !< par direct beam fraction (dimensionless)
5111    REAL(wp) ::  ratio                       !< ratio measured to potential solar radiation (dimensionless)
5112    REAL(wp) ::  rdm                         !< potential direct beam near-infrared radiation (W m-2); "potential" means clear-sky
5113    REAL(wp) ::  rdn                         !< potential diffuse near-infrared radiation (W m-2)
5114    REAL(wp) ::  rdu                         !< visible (par) direct beam radiation (W m-2)
5115    REAL(wp) ::  rdv                         !< potential visible (par) diffuse radiation (W m-2)
5116    REAL(wp) ::  rn                          !< near-infrared radiation (W m-2)
5117    REAL(wp) ::  rv                          !< visible radiation (W m-2)
5118    REAL(wp) ::  ww                          !< water absorption in the near infrared for 10 mm of precipitable water
5119
5120!
5121!-- Calculate visible (PAR) direct beam radiation
5122!-- 600 W m-2 represents average amount of par (400-700 nm wavelength)
5123!-- at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2)
5124    rdu = 600.0_wp* exp( -0.185_wp * ( pres / pres_0 ) / sinphi ) * sinphi
5125!
5126!-- Calculate potential visible diffuse radiation
5127    rdv = 0.4_wp * ( 600.0_wp - rdu ) * sinphi
5128!
5129!-- Calculate the water absorption in the-near infrared
5130    ww = 1320 * 10**( -1.195_wp + 0.4459_wp * log10( 1.0_wp / sinphi ) - 0.0345_wp * ( log10( 1.0_wp / sinphi ) )**2 )
5131!
5132!-- Calculate potential direct beam near-infrared radiation
5133    rdm = (720.0_wp * exp(-0.06_wp * (pres / pres_0) / sinphi ) - ww ) * sinphi     !< 720 = solar constant - 600
5134!
5135!-- Calculate potential diffuse near-infrared radiation
5136    rdn = 0.6_wp * ( 720 - rdm - ww ) * sinphi
5137!
5138!-- Compute visible and near-infrared radiation
5139    rv = MAX( 0.1_wp, rdu + rdv )
5140    rn = MAX( 0.01_wp, rdm + rdn )
5141!
5142!-- Compute ratio between input global radiation (here defined as solar radiation, dirict+diffuse)
5143!-- and total radiation computed here
5144    ratio = MIN( 0.89_wp, solar_rad / ( rv + rn ) )
5145!
5146!-- Calculate total visible radiation
5147    sv = ratio * rv
5148!
5149!-- Calculate fraction of par in the direct beam
5150    fv = MIN( 0.99_wp, ( 0.9_wp - ratio ) / 0.7_wp )              !< help variable
5151    fv = MAX( 0.01_wp, rdu / rv * ( 1.0_wp - fv**0.6667_wp ) )    !< fraction of par in the direct beam
5152!
5153!-- Compute direct and diffuse parts of par
5154    par_dir = fv * sv
5155    par_diff = sv - par_dir
5156
5157 END SUBROUTINE par_dir_diff
5158
5159 
5160 !-------------------------------------------------------------------
5161 !> rc_get_vpd: get vapour pressure deficit (kPa)
5162 !-------------------------------------------------------------------
5163 SUBROUTINE rc_get_vpd( temp, rh, vpd )
5164
5165!
5166!-- Input/output variables:
5167    REAL(wp), INTENT(IN) ::  temp    !< temperature (C)
5168    REAL(wp), INTENT(IN) ::  rh    !< relative humidity (%)
5169
5170    REAL(wp), INTENT(OUT) ::  vpd    !< vapour pressure deficit (kPa)
5171!
5172!-- Local variables:
5173    REAL(wp) ::  esat
5174!
5175!-- fit parameters:
5176    REAL(wp), PARAMETER ::  a1 = 6.113718e-01
5177    REAL(wp), PARAMETER ::  a2 = 4.43839e-02
5178    REAL(wp), PARAMETER ::  a3 = 1.39817e-03
5179    REAL(wp), PARAMETER ::  a4 = 2.9295e-05
5180    REAL(wp), PARAMETER ::  a5 = 2.16e-07
5181    REAL(wp), PARAMETER ::  a6 = 3.0e-09
5182!
5183!-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973)
5184    esat = a1 + a2 * temp + a3 * temp**2 + a4 * temp**3 + a5 * temp**4 + a6 * temp**5
5185    vpd  = esat * ( 1 - rh / 100 )
5186
5187 END SUBROUTINE rc_get_vpd
5188
5189
5190 !-------------------------------------------------------------------
5191 !> rc_gsoil_eff: compute effective soil conductance
5192 !-------------------------------------------------------------------
5193 SUBROUTINE rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff )
5194
5195!
5196!-- Input/output variables:
5197    INTEGER(iwp), INTENT(IN) ::  icmp          !< component index
5198    INTEGER(iwp), INTENT(IN) ::  lu            !< land use type, lu = 1,..., nlu
5199    INTEGER(iwp), INTENT(IN) ::  nwet          !< index for wetness
5200                                               !< nwet = 0 -> dry; nwet = 1 -> wet; nwet = 9 -> snow
5201                                               !< N.B. this routine cannot be called with nwet = 9,
5202                                               !< nwet = 9 should be handled outside this routine.
5203    REAL(wp), INTENT(IN) ::  sai               !< surface area index
5204    REAL(wp), INTENT(IN) ::  ust               !< friction velocity (m/s)
5205    REAL(wp), INTENT(IN) ::  t                 !< temperature (C)
5206    REAL(wp), INTENT(OUT) ::  gsoil_eff        !< effective soil conductance (m/s)
5207!
5208!-- local variables:
5209    REAL(wp) ::  rinc                          !< in canopy resistance  (s/m)
5210    REAL(wp) ::  rsoil_eff                     !< effective soil resistance (s/m)
5211!
5212!-- Soil resistance (numbers matched with lu_classes and component numbers)
5213    !     grs    ara    crp    cnf    dec    wat    urb   oth    des    ice    sav    trf    wai    med    sem
5214    REAL(wp), PARAMETER ::  rsoil(nlu_dep,ncmp) = reshape( (/ &
5215         1000.,  200.,  200.,  200.,  200., 2000.,  400., 1000., 2000., 2000., 1000.,  200., 2000.,  200.,  400., &    !< O3
5216         1000., 1000., 1000., 1000., 1000.,   10., 1000., 1000., 1000.,  500., 1000., 1000.,   10., 1000., 1000., &    !< SO2
5217         1000., 1000., 1000., 1000., 1000., 2000., 1000., 1000., 1000., 2000., 1000., 1000., 2000., 1000., 1000., &    !< NO2
5218         -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., &    !< NO
5219         100.,  100.,  100.,  100.,  100.,   10.,  100.,  100.,  100., 1000.,  100.,  100.,   10.,  100.,  100.,  &    !< NH3
5220         -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., &    !< CO
5221         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., &    !< NO3
5222         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., &    !< HNO3
5223         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., &    !< N2O5
5224         -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999., -999. /),&  !< H2O2   
5225         (/nlu_dep,ncmp/) )
5226!
5227!-- For                                          o3    so2   no2     no    nh3     co     no3    hno3   n2o5   h2o2
5228    REAL(wp), PARAMETER ::  rsoil_wet(ncmp)    = (/2000., 10. , 2000., -999., 10.  , -999., -999., -999., -999., -999./)
5229    REAL(wp), PARAMETER ::  rsoil_frozen(ncmp) = (/2000., 500., 2000., -999., 1000., -999., -999., -999., -999., -999./)
5230!
5231!-- Compute in canopy (in crop) resistance:
5232    CALL rc_rinc( lu, sai, ust, rinc )
5233!
5234!-- Check for missing deposition path:
5235    IF ( missing(rinc) )  THEN
5236       rsoil_eff = -9999.0_wp
5237    ELSE
5238!
5239!--    Frozen soil (temperature below 0):
5240       IF ( t < 0.0_wp )  THEN
5241          IF ( missing( rsoil_frozen( icmp ) ) )  THEN
5242             rsoil_eff = -9999.0_wp
5243          ELSE
5244             rsoil_eff = rsoil_frozen( icmp ) + rinc
5245          ENDIF
5246       ELSE
5247!
5248!--       Non-frozen soil; dry:
5249          IF ( nwet == 0 )  THEN
5250             IF ( missing( rsoil( lu, icmp ) ) )  THEN
5251                rsoil_eff = -9999.0_wp
5252             ELSE
5253                rsoil_eff = rsoil( lu, icmp ) + rinc
5254             ENDIF
5255!
5256!--       Non-frozen soil; wet:
5257          ELSEIF ( nwet == 1 )  THEN
5258             IF ( missing( rsoil_wet( icmp ) ) )  THEN
5259                rsoil_eff = -9999.0_wp
5260             ELSE
5261                rsoil_eff = rsoil_wet( icmp ) + rinc
5262             ENDIF
5263          ELSE
5264             message_string = 'nwet can only be 0 or 1'
5265             CALL message( 'rc_gsoil_eff', 'CM0460', 1, 2, 0, 6, 0 )
5266          ENDIF
5267       ENDIF
5268    ENDIF
5269!
5270!-- Compute conductance:
5271    IF ( rsoil_eff > 0.0_wp )  THEN
5272       gsoil_eff = 1.0_wp / rsoil_eff
5273    ELSE
5274       gsoil_eff = 0.0_wp
5275    ENDIF
5276
5277 END SUBROUTINE rc_gsoil_eff
5278
5279
5280 !-------------------------------------------------------------------
5281 !> rc_rinc: compute in canopy (or in crop) resistance
5282 !> van Pul and Jacobs, 1993, BLM
5283 !-------------------------------------------------------------------
5284 SUBROUTINE rc_rinc( lu, sai, ust, rinc )
5285
5286!
5287!-- Input/output variables:
5288    INTEGER(iwp), INTENT(IN) ::  lu          !< land use class, lu = 1, ..., nlu
5289
5290    REAL(wp), INTENT(IN) ::  sai             !< surface area index
5291    REAL(wp), INTENT(IN) ::  ust             !< friction velocity (m/s)
5292
5293    REAL(wp), INTENT(OUT) ::  rinc           !< in canopy resistance (s/m)
5294!
5295!-- b = empirical constant for computation of rinc (in canopy resistance) (= 14 m-1 or -999 if not applicable)
5296!-- h = vegetation height (m)                     gra  ara crop con dec wat   urb   oth   des   ice   sav   trf  wai  med semi
5297    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  b = (/ -999, 14, 14, 14, 14, -999, -999, -999, -999, -999, -999, 14, -999,  &
5298         14, 14 /)
5299    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  h = (/ -999, 1,  1,  20, 20, -999, -999, -999, -999, -999, -999, 20, -999,  &
5300         1 ,  1 /)
5301!
5302!-- Compute Rinc only for arable land, perm. crops, forest; otherwise Rinc = 0:
5303    IF ( b(lu) > 0.0_wp )  THEN
5304!       !
5305!--    Check for u* > 0 (otherwise denominator = 0):
5306       IF ( ust > 0.0_wp )  THEN
5307          rinc = b(lu) * h(lu) * sai/ust
5308       ELSE
5309          rinc = 1000.0_wp
5310       ENDIF
5311    ELSE
5312       IF ( lu == ilu_grass .OR. lu == ilu_other  )  THEN
5313          rinc = -999.0_wp     !< no deposition path for grass, other, and semi-natural
5314       ELSE
5315          rinc = 0.0_wp        !< no in-canopy resistance
5316       ENDIF
5317    ENDIF
5318
5319 END SUBROUTINE rc_rinc
5320
5321
5322 !-------------------------------------------------------------------
5323 !> rc_rctot: compute total canopy (or surface) resistance Rc
5324 !-------------------------------------------------------------------
5325 SUBROUTINE rc_rctot( gstom, gsoil_eff, gw, gc_tot, rc_tot )
5326
5327!
5328!-- Input/output variables:
5329    REAL(wp), INTENT(IN) ::  gstom         !< stomatal conductance (s/m)
5330    REAL(wp), INTENT(IN) ::  gsoil_eff     !< effective soil conductance (s/m)
5331    REAL(wp), INTENT(IN) ::  gw            !< external leaf conductance (s/m)
5332
5333    REAL(wp), INTENT(OUT) ::  gc_tot       !< total canopy conductance (m/s)
5334    REAL(wp), INTENT(OUT) ::  rc_tot       !< total canopy resistance Rc (s/m)
5335!
5336!-- Total conductance:
5337    gc_tot = gstom + gsoil_eff + gw
5338!
5339!-- Total resistance (note: gw can be negative, but no total emission allowed here):
5340    IF ( gc_tot <= 0.0_wp .OR. gw < 0.0_wp )  THEN
5341       rc_tot = -9999.0_wp
5342    ELSE
5343       rc_tot = 1.0_wp / gc_tot
5344    ENDIF
5345
5346 END SUBROUTINE rc_rctot
5347
5348
5349 !-------------------------------------------------------------------
5350 !> rc_comp_point_rc_eff: calculate the effective resistance Rc
5351 !> based on one or more compensation points
5352 !-------------------------------------------------------------------
5353 !> NH3rc (see depac v3.6 is based on Avero workshop Marc Sutton. p. 173.
5354 !> Sutton 1998 AE 473-480)
5355 !>
5356 !> Documentation by Ferd Sauter, 2008; see also documentation block in header of depac subroutine.
5357 !> FS 2009-01-29: variable names made consistent with DEPAC
5358 !> FS 2009-03-04: use total compensation point
5359 !>
5360 !> C: with total compensation point   ! D: approximation of C
5361 !>                                    !    with classical approach
5362 !>  zr --------- Catm                 !  zr --------- Catm   
5363 !>         |                          !         |       
5364 !>         Ra                         !         Ra     
5365 !>         |                          !         |       
5366 !>         Rb                         !         Rb     
5367 !>         |                          !         |       
5368 !>  z0 --------- Cc                   !  z0 --------- Cc
5369 !>         |                          !         |             
5370 !>        Rc                          !        Rc_eff         
5371 !>         |                          !         |             
5372 !>     --------- Ccomp_tot            !     --------- C=0
5373 !>
5374 !>
5375 !> The effective Rc is defined such that instead of using
5376 !>
5377 !>   F = -vd*[Catm - Ccomp_tot]                                    (1)
5378 !>
5379 !> we can use the 'normal' flux formula
5380 !>
5381 !>   F = -vd'*Catm,                                                (2)
5382 !>
5383 !> with vd' = 1/(Ra + Rb + Rc')                                    (3)
5384 !>
5385 !> and Rc' the effective Rc (rc_eff).
5386 !>                                                (Catm - Ccomp_tot)
5387 !> vd'*Catm = vd*(Catm - Ccomp_tot) <=> vd' = vd* ------------------
5388 !>                                                      Catm
5389 !>
5390 !>                                        (Catm - Ccomp_tot)
5391 !> 1/(Ra + Rb + Rc') = (1/Ra + Rb + Rc) * ------------------
5392 !>                                              Catm
5393 !>
5394 !>                                          Catm
5395 !> (Ra + Rb + Rc') = (Ra + Rb + Rc) * ------------------
5396 !>                                     (Catm - Ccomp_tot)
5397 !>
5398 !>                              Catm
5399 !> Rc' = (Ra + Rb + Rc) * ------------------ - Ra - Rb
5400 !>                        (Catm - Ccomp_tot)
5401 !>
5402 !>                        Catm                           Catm
5403 !> Rc' = (Ra + Rb) [------------------ - 1 ] + Rc * ------------------
5404 !>                  (Catm - Ccomp_tot)              (Catm - Ccomp_tot)
5405 !>
5406 !> Rc' = [(Ra + Rb)*Ccomp_tot + Rc*Catm ] / (Catm - Ccomp_tot)
5407 !>
5408 ! -------------------------------------------------------------------------------------------
5409! SUBROUTINE rc_comp_point_rc_eff( ccomp_tot, conc_ijk_ugm3, ra, rb, rc_tot, rc_eff )
5410!
5411!
5412!!-- Input/output variables:
5413!    REAL(wp), INTENT(IN) ::  ccomp_tot     !< total compensation point (weighed average of separate compensation points) (ug/m3)
5414!    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3 !< atmospheric concentration (ug/m3) above Catm
5415!    REAL(wp), INTENT(IN) ::  ra            !< aerodynamic resistance (s/m)
5416!    REAL(wp), INTENT(IN) ::  rb            !< boundary layer resistance (s/m)
5417!    REAL(wp), INTENT(IN) ::  rc_tot        !< total canopy resistance (s/m)
5418!
5419!    REAL(wp), INTENT(OUT) ::  rc_eff       !< effective total canopy resistance (s/m)
5420!
5421!    !
5422!!-- Compute effective resistance:
5423!    IF (  ccomp_tot == 0.0_wp )  THEN
5424!       !
5425!!--    trace with no compensiation point ( or compensation point equal to zero)
5426!       rc_eff = rc_tot
5427!
5428!    ELSE IF ( ccomp_tot > 0.0_wp .AND. ( abs( conc_ijk_ugm3 - ccomp_tot ) < 1.e-8 ) )  THEN
5429!       !
5430!!--   surface concentration (almost) equal to atmospheric concentration
5431!!--    no exchange between surface and atmosphere, infinite RC --> vd=0
5432!       rc_eff = 9999999999.0_wp
5433!
5434!    ELSE IF ( ccomp_tot > 0.0_wp )  THEN
5435!       !
5436!!--    compensation point available, calculate effective resistance
5437!       rc_eff = ( ( ra + rb ) * ccomp_tot + rc_tot * conc_ijk_ugm3 ) / ( conc_ijk_ugm3 - ccomp_tot )
5438!
5439!    ELSE
5440!       rc_eff = -999.0_wp
5441!       message_string = 'This should not be possible, check ccomp_tot'
5442!       CALL message( 'rc_comp_point_rc_eff', 'CM0461', 1, 2, 0, 6, 0 )
5443!    ENDIF
5444!
5445!    RETURN
5446!   
5447! END SUBROUTINE rc_comp_point_rc_eff
5448
5449
5450 !-------------------------------------------------------------------
5451 !> missing: check for data that correspond with a missing deposition path
5452 !>          this data is represented by -999
5453 !-------------------------------------------------------------------
5454 LOGICAL function missing( x )
5455
5456    REAL(wp), INTENT(IN) ::  x
5457
5458!
5459!-- bandwidth for checking (in)equalities of floats
5460    REAL(wp), PARAMETER :: eps = 1.0e-5
5461
5462    missing = (abs(x + 999.0_wp) <= eps)
5463
5464 END function missing
5465