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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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