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

Last change on this file since 3719 was 3719, checked in by kanani, 4 years ago

Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)

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