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

Last change on this file since 3880 was 3880, checked in by knoop, 2 years ago

Moved chem_prognostic_equations into module_interface

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