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

Last change on this file since 4346 was 4346, checked in by motisi, 23 months ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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