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

Last change on this file since 4238 was 4230, checked in by suehring, 5 years ago

Several bugfixes: profile initialization of chemical species in restart runs; Runge-Kutta tendency array not initialized in chemistry model in restart runs; fixed determination of time indices for chemical emissions (introduced with commit -4227); Update chemistry profiles in offline nesting; initialize canopy resistances for greened building walls (even if green fraction is zero)

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