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

Last change on this file since 4273 was 4273, checked in by monakurppa, 18 months ago

Add logical switched nesting_chem and nesting_offline_chem

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