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

Last change on this file since 4442 was 4442, checked in by suehring, 4 years ago

last commit documented

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