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

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

Implementation of a monotonic flux limiter for vertical advection term in Wicker-Skamarock scheme. The flux limiter is currently only applied for passive scalars (passive scalar, chemical species, aerosols) within the region up to the highest topography, in order to avoid the built-up of large concentrations within poorly resolved cavities in urban environments. To enable the limiter monotonic_limiter_z = .T. must be set. Note, the limiter is currently only implemented for the cache-optimized version of advec_ws. Further changes in offline nesting: Set boundary condition for w at nzt+1 at all lateral boundaries (even though these won't enter the numerical solution), in order to avoid high vertical velocities in the run-control file which might built-up due to the mass-conservation; bugfix in offline nesting for chemical species

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