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

Last change on this file since 3877 was 3877, checked in by knoop, 2 years ago

Added chem_actions module interface to chemistry_model_mod and moved photolysis_control call into it

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