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

Last change on this file since 4550 was 4550, checked in by raasch, 11 months ago

bugfix for reading local restart data

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