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

Last change on this file since 3848 was 3848, checked in by forkel, 2 years ago

some formatting

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