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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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