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

Last change on this file since 3930 was 3930, checked in by forkel, 3 years ago

changed subroutine name from chem_non_transport_physics to chem_non_advective_processes

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