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

Last change on this file since 4457 was 4457, checked in by raasch, 4 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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