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

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

Restore accidantly removed limitation to positive values

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