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

Last change on this file since 4268 was 4268, checked in by schwenkel, 21 months ago

Introducing module interface for boundary conditions and move module specific boundary conditions into their modules

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