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

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

Added chem_non_transport_physics module interface to chemistry_model_mod and moved respective calls into it

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