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

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

2D output of emission fluxes in chemistry; revise check for multigrid coarsening and give a warning instead of an error

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