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

Last change on this file since 4227 was 4227, checked in by gronemeier, 22 months ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

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