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

Last change on this file since 3898 was 3898, checked in by suehring, 2 years ago

remove un-used variable

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