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

Last change on this file since 3685 was 3685, checked in by knoop, 3 years ago

Some interface calls moved to module_interface + cleanup

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