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

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

Updates from chemistriy branched merged into trunk: code cleaning and formatting, code structure optimizations

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