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

Last change on this file since 3929 was 3929, checked in by banzhafs, 3 years ago

Correct/complete module_interface introduction for chemistry model and bug fix in chem_depo subroutine

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