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

Last change on this file since 4372 was 4372, checked in by banzhafs, 19 months ago

added namelist flag emiss_read_legacy_mode to allow concurrent

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