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

Last change on this file since 3862 was 3862, checked in by banzhafs, 2 years ago

some formatting in deposition code of chemistry_model_mod

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