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

Last change on this file since 4292 was 4292, checked in by banzhafs, 17 months ago

Bugfix for r4290

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