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

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

Bugfix: replace simulated_time by time_since_reference_point where required

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