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

Last change on this file since 3687 was 3687, checked in by knoop, 5 years ago

Moved all user routunes that are dependencies of the PALM core only, to user_module.f90
The files that formerly contained these routines, have been deleted.
Also module_interface routines for init_mask and last_actions have been added.

  • 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 3687 2019-01-22 10:42:06Z 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    IMPLICIT NONE
2633
2634    CHARACTER (LEN=*) ::  mode   !<
2635
2636    INTEGER(iwp) ::  i    !< running index on x-axis
2637    INTEGER(iwp) ::  j    !< running index on y-axis
2638    INTEGER(iwp) ::  k    !< vertical index counter
2639    INTEGER(iwp) ::  sr   !< statistical region
2640    INTEGER(iwp) ::  tn   !< thread number
2641    INTEGER(iwp) ::  n    !<
2642    INTEGER(iwp) ::  m    !<
2643    INTEGER(iwp) ::  lpr  !< running index chem spcs
2644
2645    !    REAL(wp),                                                                                      &
2646    !    DIMENSION(dots_num_palm+1:dots_max) ::                                                         &
2647    !          ts_value_l   !<
2648
2649    IF ( mode == 'profiles' )  THEN
2650       !
2651       !--    Sample on how to calculate horizontally averaged profiles of user-
2652       !--    defined quantities. Each quantity is identified by the index
2653       !--    "pr_palm+#" where "#" is an integer starting from 1. These
2654       !--    user-profile-numbers must also be assigned to the respective strings
2655       !--    given by data_output_pr_cs in routine user_check_data_output_pr.
2656       !--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
2657       !                       w*pt*), dim-4 = statistical region.
2658
2659       !$OMP DO
2660       DO  i = nxl, nxr
2661          DO  j = nys, nyn
2662             DO  k = nzb, nzt+1
2663                DO lpr = 1, cs_pr_count
2664
2665                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
2666                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
2667                        rmask(j,i,sr)  *                                   &
2668                        MERGE( 1.0_wp, 0.0_wp,                             &
2669                        BTEST( wall_flags_0(k,j,i), 22 ) )
2670                ENDDO
2671             ENDDO
2672          ENDDO
2673       ENDDO
2674    ELSEIF ( mode == 'time_series' )  THEN
2675!      @todo
2676    ENDIF
2677
2678 END SUBROUTINE chem_statistics
2679
2680 !
2681 !------------------------------------------------------------------------------!
2682 !
2683 ! Description:
2684 ! ------------
2685 !> Subroutine for swapping of timelevels for chemical species
2686 !> called out from subroutine swap_timelevel
2687 !------------------------------------------------------------------------------!
2688
2689 SUBROUTINE chem_swap_timelevel( level )
2690
2691    IMPLICIT NONE
2692
2693    INTEGER(iwp), INTENT(IN) ::  level
2694    !
2695    !-- local variables
2696    INTEGER(iwp)             ::  lsp
2697
2698
2699    IF ( level == 0 ) THEN
2700       DO lsp=1, nvar                                       
2701          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
2702          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
2703       ENDDO
2704    ELSE
2705       DO lsp=1, nvar                                       
2706          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
2707          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
2708       ENDDO
2709    ENDIF
2710
2711    RETURN
2712 END SUBROUTINE chem_swap_timelevel
2713 !
2714 !------------------------------------------------------------------------------!
2715 !
2716 ! Description:
2717 ! ------------
2718 !> Subroutine to write restart data for chemistry model
2719 !------------------------------------------------------------------------------!
2720 SUBROUTINE chem_wrd_local
2721
2722    IMPLICIT NONE
2723
2724    INTEGER(iwp) ::  lsp  !<
2725
2726    DO  lsp = 1, nspec
2727       CALL wrd_write_string( TRIM( chem_species(lsp)%name ) )
2728       WRITE ( 14 )  chem_species(lsp)%conc
2729       CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
2730       WRITE ( 14 )  chem_species(lsp)%conc_av
2731    ENDDO
2732
2733 END SUBROUTINE chem_wrd_local
2734
2735
2736
2737 !-------------------------------------------------------------------------------!
2738 !
2739 ! Description:
2740 ! ------------
2741 !> Subroutine to calculate the deposition of gases and PMs. For now deposition
2742 !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
2743 !> considered. The deposition of particles is derived following Zhang et al.,
2744 !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
2745 !>     
2746 !> @TODO: Consider deposition on vertical surfaces   
2747 !> @TODO: Consider overlaying horizontal surfaces
2748 !> @TODO: Consider resolved vegetation   
2749 !> @TODO: Check error messages
2750 !-------------------------------------------------------------------------------!
2751
2752 SUBROUTINE chem_depo( i, j )
2753
2754    USE control_parameters,                                                 &   
2755         ONLY:  dt_3d, intermediate_timestep_count, latitude
2756
2757    USE arrays_3d,                                                          &
2758         ONLY:  dzw, rho_air_zw
2759
2760    USE date_and_time_mod,                                                  &
2761         ONLY:  day_of_year
2762
2763    USE surface_mod,                                                        &
2764         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h,        &
2765         surf_type, surf_usm_h
2766
2767    USE radiation_model_mod,                                                &
2768         ONLY:  zenith
2769
2770
2771
2772    IMPLICIT NONE
2773
2774    INTEGER(iwp), INTENT(IN) ::  i
2775    INTEGER(iwp), INTENT(IN) ::  j
2776
2777    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
2778    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
2779    INTEGER(iwp) ::  lu                            !< running index for landuse classes
2780    INTEGER(iwp) ::  lu_palm                       !< index of PALM LSM vegetation_type at current surface element
2781    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
2782    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current surface element
2783    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current surface element
2784    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current surface element
2785    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface element
2786    INTEGER(iwp) ::  lu_dep                        !< matching DEPAC LU to lu_palm
2787    INTEGER(iwp) ::  lup_dep                       !< matching DEPAC LU to lup_palm
2788    INTEGER(iwp) ::  luw_dep                       !< matching DEPAC LU to luw_palm
2789    INTEGER(iwp) ::  luu_dep                       !< matching DEPAC LU to luu_palm
2790    INTEGER(iwp) ::  lug_dep                       !< matching DEPAC LU to lug_palm
2791    INTEGER(iwp) ::  lud_dep                       !< matching DEPAC LU to lud_palm
2792    INTEGER(iwp) ::  m                             !< index for horizontal surfaces
2793
2794    INTEGER(iwp) ::  pspec                         !< running index
2795    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
2796
2797    !
2798    !-- Vegetation                                              !<Match to DEPAC
2799    INTEGER(iwp) ::  ind_lu_user = 0                         !< --> ERROR 
2800    INTEGER(iwp) ::  ind_lu_b_soil = 1                       !< --> ilu_desert
2801    INTEGER(iwp) ::  ind_lu_mixed_crops = 2                  !< --> ilu_arable
2802    INTEGER(iwp) ::  ind_lu_s_grass = 3                      !< --> ilu_grass
2803    INTEGER(iwp) ::  ind_lu_ev_needle_trees = 4              !< --> ilu_coniferous_forest
2804    INTEGER(iwp) ::  ind_lu_de_needle_trees = 5              !< --> ilu_coniferous_forest
2805    INTEGER(iwp) ::  ind_lu_ev_broad_trees = 6               !< --> ilu_tropical_forest
2806    INTEGER(iwp) ::  ind_lu_de_broad_trees = 7               !< --> ilu_deciduous_forest
2807    INTEGER(iwp) ::  ind_lu_t_grass = 8                      !< --> ilu_grass
2808    INTEGER(iwp) ::  ind_lu_desert = 9                       !< --> ilu_desert
2809    INTEGER(iwp) ::  ind_lu_tundra = 10                      !< --> ilu_other
2810    INTEGER(iwp) ::  ind_lu_irr_crops = 11                   !< --> ilu_arable
2811    INTEGER(iwp) ::  ind_lu_semidesert = 12                  !< --> ilu_other
2812    INTEGER(iwp) ::  ind_lu_ice = 13                         !< --> ilu_ice
2813    INTEGER(iwp) ::  ind_lu_marsh = 14                       !< --> ilu_other
2814    INTEGER(iwp) ::  ind_lu_ev_shrubs = 15                   !< --> ilu_mediterrean_scrub
2815    INTEGER(iwp) ::  ind_lu_de_shrubs = 16                   !< --> ilu_mediterrean_scrub
2816    INTEGER(iwp) ::  ind_lu_mixed_forest = 17                !< --> ilu_coniferous_forest (ave(decid+conif))
2817    INTEGER(iwp) ::  ind_lu_intrup_forest = 18               !< --> ilu_other (ave(other+decid))
2818
2819    !
2820    !-- Water
2821    INTEGER(iwp) ::  ind_luw_user = 0                        !< --> ERROR
2822    INTEGER(iwp) ::  ind_luw_lake = 1                        !< --> ilu_water_inland
2823    INTEGER(iwp) ::  ind_luw_river = 2                       !< --> ilu_water_inland
2824    INTEGER(iwp) ::  ind_luw_ocean = 3                       !< --> ilu_water_sea
2825    INTEGER(iwp) ::  ind_luw_pond = 4                        !< --> ilu_water_inland
2826    INTEGER(iwp) ::  ind_luw_fountain = 5                    !< --> ilu_water_inland
2827
2828    !
2829    !-- Pavement
2830    INTEGER(iwp) ::  ind_lup_user = 0                        !< --> ERROR
2831    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !< --> ilu_desert
2832    INTEGER(iwp) ::  ind_lup_asph = 2                        !< --> ilu_desert
2833    INTEGER(iwp) ::  ind_lup_conc = 3                        !< --> ilu_desert
2834    INTEGER(iwp) ::  ind_lup_sett = 4                        !< --> ilu_desert
2835    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !< --> ilu_desert
2836    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !< --> ilu_desert
2837    INTEGER(iwp) ::  ind_lup_metal = 7                       !< --> ilu_desert
2838    INTEGER(iwp) ::  ind_lup_wood = 8                        !< --> ilu_desert
2839    INTEGER(iwp) ::  ind_lup_gravel = 9                      !< --> ilu_desert
2840    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !< --> ilu_desert
2841    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !< --> ilu_desert
2842    INTEGER(iwp) ::  ind_lup_woodchips = 12                  !< --> ilu_desert
2843    INTEGER(iwp) ::  ind_lup_tartan = 13                     !< --> ilu_desert
2844    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !< --> ilu_desert
2845    INTEGER(iwp) ::  ind_lup_clay = 15                       !< --> ilu_desert
2846
2847
2848    !
2849    !-- Particle parameters according to the respective aerosol classes (PM25, PM10)
2850    INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
2851    INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
2852    INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
2853
2854    INTEGER(iwp) ::  part_type
2855
2856    INTEGER(iwp) ::  nwet           !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
2857
2858    REAL(wp) ::  dt_chem
2859    REAL(wp) ::  dh               
2860    REAL(wp) ::  inv_dh
2861    REAL(wp) ::  dt_dh
2862
2863    REAL(wp) ::  dens              !< density at layer k at i,j 
2864    REAL(wp) ::  r_aero_lu         !< aerodynamic resistance (s/m) at current surface element
2865    REAL(wp) ::  ustar_lu          !< ustar at current surface element
2866    REAL(wp) ::  z0h_lu            !< roughness length for heat at current surface element
2867    REAL(wp) ::  glrad             !< rad_sw_in at current surface element
2868    REAL(wp) ::  ppm_to_ugm3       !< conversion factor
2869    REAL(wp) ::  relh              !< relative humidity at current surface element
2870    REAL(wp) ::  lai               !< leaf area index at current surface element
2871    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
2872
2873    REAL(wp) ::  slinnfac       
2874    REAL(wp) ::  visc              !< Viscosity
2875    REAL(wp) ::  vs                !< Sedimentation velocity
2876    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
2877    REAL(wp) ::  Rs                !< Sedimentaion resistance (s/m)
2878    REAL(wp) ::  Rb                !< quasi-laminar boundary layer resistance (s/m)
2879    REAL(wp) ::  Rc_tot            !< total canopy resistance Rc (s/m)
2880
2881    REAL(wp) ::  catm              !< gasconc. at i, j, k in ug/m3
2882    REAL(wp) ::  diffc             !< diffusivity
2883
2884
2885    REAL(wp), DIMENSION(nspec) ::  bud_now_lu       !< budget for LSM vegetation type at current surface element
2886    REAL(wp), DIMENSION(nspec) ::  bud_now_lup      !< budget for LSM pavement type at current surface element
2887    REAL(wp), DIMENSION(nspec) ::  bud_now_luw      !< budget for LSM water type at current surface element
2888    REAL(wp), DIMENSION(nspec) ::  bud_now_luu      !< budget for USM walls/roofs at current surface element
2889    REAL(wp), DIMENSION(nspec) ::  bud_now_lug      !< budget for USM green surfaces at current surface element
2890    REAL(wp), DIMENSION(nspec) ::  bud_now_lud      !< budget for USM windows at current surface element
2891    REAL(wp), DIMENSION(nspec) ::  bud_now          !< overall budget at current surface element
2892    REAL(wp), DIMENSION(nspec) ::  cc               !< concentration i,j,k
2893    REAL(wp), DIMENSION(nspec) ::  ccomp_tot        !< total compensation point (ug/m3)
2894    !< For now kept to zero for all species!
2895
2896    REAL(wp) ::  ttemp          !< temperatur at i,j,k
2897    REAL(wp) ::  ts             !< surface temperatur in degrees celsius
2898    REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
2899
2900    !
2901    !-- Particle parameters (PM10 (1), PM25 (2))
2902    !-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
2903    !-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
2904    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/ &
2905         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &  !<  1
2906         0.7e-6_wp, 1.14e3_wp, 1.082_wp &   !<  2
2907         /), (/ 3, 2 /) )
2908
2909
2910    LOGICAL ::  match_lsm     !< flag indicating natural-type surface
2911    LOGICAL ::  match_usm     !< flag indicating urban-type surface
2912
2913
2914    !
2915    !-- List of names of possible tracers
2916    CHARACTER(len=*), PARAMETER ::  pspecnames(nposp) = (/ &
2917         'NO2           ', &    !< NO2
2918         'NO            ', &    !< NO
2919         'O3            ', &    !< O3
2920         'CO            ', &    !< CO
2921         'form          ', &    !< FORM
2922         'ald           ', &    !< ALD
2923         'pan           ', &    !< PAN
2924         'mgly          ', &    !< MGLY
2925         'par           ', &    !< PAR
2926         'ole           ', &    !< OLE
2927         'eth           ', &    !< ETH
2928         'tol           ', &    !< TOL
2929         'cres          ', &    !< CRES
2930         'xyl           ', &    !< XYL
2931         'SO4a_f        ', &    !< SO4a_f
2932         'SO2           ', &    !< SO2
2933         'HNO2          ', &    !< HNO2
2934         'CH4           ', &    !< CH4
2935         'NH3           ', &    !< NH3
2936         'NO3           ', &    !< NO3
2937         'OH            ', &    !< OH
2938         'HO2           ', &    !< HO2
2939         'N2O5          ', &    !< N2O5
2940         'SO4a_c        ', &    !< SO4a_c
2941         'NH4a_f        ', &    !< NH4a_f
2942         'NO3a_f        ', &    !< NO3a_f
2943         'NO3a_c        ', &    !< NO3a_c
2944         'C2O3          ', &    !< C2O3
2945         'XO2           ', &    !< XO2
2946         'XO2N          ', &    !< XO2N
2947         'cro           ', &    !< CRO
2948         'HNO3          ', &    !< HNO3
2949         'H2O2          ', &    !< H2O2
2950         'iso           ', &    !< ISO
2951         'ispd          ', &    !< ISPD
2952         'to2           ', &    !< TO2
2953         'open          ', &    !< OPEN
2954         'terp          ', &    !< TERP
2955         'ec_f          ', &    !< EC_f
2956         'ec_c          ', &    !< EC_c
2957         'pom_f         ', &    !< POM_f
2958         'pom_c         ', &    !< POM_c
2959         'ppm_f         ', &    !< PPM_f
2960         'ppm_c         ', &    !< PPM_c
2961         'na_ff         ', &    !< Na_ff
2962         'na_f          ', &    !< Na_f
2963         'na_c          ', &    !< Na_c
2964         'na_cc         ', &    !< Na_cc
2965         'na_ccc        ', &    !< Na_ccc
2966         'dust_ff       ', &    !< dust_ff
2967         'dust_f        ', &    !< dust_f
2968         'dust_c        ', &    !< dust_c
2969         'dust_cc       ', &    !< dust_cc
2970         'dust_ccc      ', &    !< dust_ccc
2971         'tpm10         ', &    !< tpm10
2972         'tpm25         ', &    !< tpm25
2973         'tss           ', &    !< tss
2974         'tdust         ', &    !< tdust
2975         'tc            ', &    !< tc
2976         'tcg           ', &    !< tcg
2977         'tsoa          ', &    !< tsoa
2978         'tnmvoc        ', &    !< tnmvoc
2979         'SOxa          ', &    !< SOxa
2980         'NOya          ', &    !< NOya
2981         'NHxa          ', &    !< NHxa
2982         'NO2_obs       ', &    !< NO2_obs
2983         'tpm10_biascorr', &    !< tpm10_biascorr
2984         'tpm25_biascorr', &    !< tpm25_biascorr
2985         'O3_biascorr   ' /)    !< o3_biascorr
2986
2987
2988    !
2989    !-- tracer mole mass:
2990    REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
2991         xm_O * 2 + xm_N, &                         !< NO2
2992         xm_O + xm_N, &                             !< NO
2993         xm_O * 3, &                                !< O3
2994         xm_C + xm_O, &                             !< CO
2995         xm_H * 2 + xm_C + xm_O, &                  !< FORM
2996         xm_H * 3 + xm_C * 2 + xm_O, &              !< ALD
2997         xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, &   !< PAN
2998         xm_H * 4 + xm_C * 3 + xm_O * 2, &          !< MGLY
2999         xm_H * 3 + xm_C, &                         !< PAR
3000         xm_H * 3 + xm_C * 2, &                     !< OLE
3001         xm_H * 4 + xm_C * 2, &                     !< ETH
3002         xm_H * 8 + xm_C * 7, &                     !< TOL
3003         xm_H * 8 + xm_C * 7 + xm_O, &              !< CRES
3004         xm_H * 10 + xm_C * 8, &                    !< XYL
3005         xm_S + xm_O * 4, &                         !< SO4a_f
3006         xm_S + xm_O * 2, &                         !< SO2
3007         xm_H + xm_O * 2 + xm_N, &                  !< HNO2
3008         xm_H * 4 + xm_C, &                         !< CH4
3009         xm_H * 3 + xm_N, &                         !< NH3
3010         xm_O * 3 + xm_N, &                         !< NO3
3011         xm_H + xm_O, &                             !< OH
3012         xm_H + xm_O * 2, &                         !< HO2
3013         xm_O * 5 + xm_N * 2, &                     !< N2O5
3014         xm_S + xm_O * 4, &                         !< SO4a_c
3015         xm_H * 4 + xm_N, &                         !< NH4a_f
3016         xm_O * 3 + xm_N, &                         !< NO3a_f
3017         xm_O * 3 + xm_N, &                         !< NO3a_c
3018         xm_C * 2 + xm_O * 3, &                     !< C2O3
3019         xm_dummy, &                                !< XO2
3020         xm_dummy, &                                !< XO2N
3021         xm_dummy, &                                !< CRO
3022         xm_H + xm_O * 3 + xm_N, &                  !< HNO3
3023         xm_H * 2 + xm_O * 2, &                     !< H2O2
3024         xm_H * 8 + xm_C * 5, &                     !< ISO
3025         xm_dummy, &                                !< ISPD
3026         xm_dummy, &                                !< TO2
3027         xm_dummy, &                                !< OPEN
3028         xm_H * 16 + xm_C * 10, &                   !< TERP
3029         xm_dummy, &                                !< EC_f
3030         xm_dummy, &                                !< EC_c
3031         xm_dummy, &                                !< POM_f
3032         xm_dummy, &                                !< POM_c
3033         xm_dummy, &                                !< PPM_f
3034         xm_dummy, &                                !< PPM_c
3035         xm_Na, &                                   !< Na_ff
3036         xm_Na, &                                   !< Na_f
3037         xm_Na, &                                   !< Na_c
3038         xm_Na, &                                   !< Na_cc
3039         xm_Na, &                                   !< Na_ccc
3040         xm_dummy, &                                !< dust_ff
3041         xm_dummy, &                                !< dust_f
3042         xm_dummy, &                                !< dust_c
3043         xm_dummy, &                                !< dust_cc
3044         xm_dummy, &                                !< dust_ccc
3045         xm_dummy, &                                !< tpm10
3046         xm_dummy, &                                !< tpm25
3047         xm_dummy, &                                !< tss
3048         xm_dummy, &                                !< tdust
3049         xm_dummy, &                                !< tc
3050         xm_dummy, &                                !< tcg
3051         xm_dummy, &                                !< tsoa
3052         xm_dummy, &                                !< tnmvoc
3053         xm_dummy, &                                !< SOxa
3054         xm_dummy, &                                !< NOya
3055         xm_dummy, &                                !< NHxa
3056         xm_O * 2 + xm_N, &                         !< NO2_obs
3057         xm_dummy, &                                !< tpm10_biascorr
3058         xm_dummy, &                                !< tpm25_biascorr
3059         xm_O * 3 /)                                !< o3_biascorr
3060
3061
3062    !
3063    !-- Initialize surface element m
3064    m = 0
3065    k = 0
3066    !
3067    !-- LSM or USM surface present at i,j:
3068    !-- Default surfaces are NOT considered for deposition
3069    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
3070    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
3071
3072
3073    !
3074    !--For LSM surfaces
3075
3076    IF ( match_lsm )  THEN
3077
3078       !
3079       !-- Get surface element information at i,j:
3080       m = surf_lsm_h%start_index(j,i)
3081       k = surf_lsm_h%k(m)
3082
3083       !
3084       !-- Get needed variables for surface element m
3085       ustar_lu  = surf_lsm_h%us(m)
3086       z0h_lu    = surf_lsm_h%z0h(m)
3087       r_aero_lu = surf_lsm_h%r_a(m)
3088       glrad     = surf_lsm_h%rad_sw_in(m)
3089       lai = surf_lsm_h%lai(m)
3090       sai = lai + 1
3091
3092       !
3093       !-- For small grid spacing neglect R_a
3094       IF ( dzw(k) <= 1.0 )  THEN
3095          r_aero_lu = 0.0_wp
3096       ENDIF
3097
3098       !
3099       !--Initialize lu's
3100       lu_palm = 0
3101       lu_dep = 0
3102       lup_palm = 0
3103       lup_dep = 0
3104       luw_palm = 0
3105       luw_dep = 0
3106
3107       !
3108       !--Initialize budgets
3109       bud_now_lu  = 0.0_wp
3110       bud_now_lup = 0.0_wp
3111       bud_now_luw = 0.0_wp
3112
3113
3114       !
3115       !-- Get land use for i,j and assign to DEPAC lu
3116       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
3117          lu_palm = surf_lsm_h%vegetation_type(m)
3118          IF ( lu_palm == ind_lu_user )  THEN
3119             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
3120             CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
3121          ELSEIF ( lu_palm == ind_lu_b_soil )  THEN
3122             lu_dep = 9
3123          ELSEIF ( lu_palm == ind_lu_mixed_crops )  THEN
3124             lu_dep = 2
3125          ELSEIF ( lu_palm == ind_lu_s_grass )  THEN
3126             lu_dep = 1
3127          ELSEIF ( lu_palm == ind_lu_ev_needle_trees )  THEN
3128             lu_dep = 4
3129          ELSEIF ( lu_palm == ind_lu_de_needle_trees )  THEN
3130             lu_dep = 4
3131          ELSEIF ( lu_palm == ind_lu_ev_broad_trees )  THEN
3132             lu_dep = 12
3133          ELSEIF ( lu_palm == ind_lu_de_broad_trees )  THEN
3134             lu_dep = 5
3135          ELSEIF ( lu_palm == ind_lu_t_grass )  THEN
3136             lu_dep = 1
3137          ELSEIF ( lu_palm == ind_lu_desert )  THEN
3138             lu_dep = 9
3139          ELSEIF ( lu_palm == ind_lu_tundra )  THEN
3140             lu_dep = 8
3141          ELSEIF ( lu_palm == ind_lu_irr_crops )  THEN
3142             lu_dep = 2
3143          ELSEIF ( lu_palm == ind_lu_semidesert )  THEN
3144             lu_dep = 8
3145          ELSEIF ( lu_palm == ind_lu_ice )  THEN
3146             lu_dep = 10
3147          ELSEIF ( lu_palm == ind_lu_marsh )  THEN
3148             lu_dep = 8
3149          ELSEIF ( lu_palm == ind_lu_ev_shrubs )  THEN
3150             lu_dep = 14
3151          ELSEIF ( lu_palm == ind_lu_de_shrubs )  THEN
3152             lu_dep = 14
3153          ELSEIF ( lu_palm == ind_lu_mixed_forest )  THEN
3154             lu_dep = 4
3155          ELSEIF ( lu_palm == ind_lu_intrup_forest )  THEN
3156             lu_dep = 8     
3157          ENDIF
3158       ENDIF
3159
3160       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
3161          lup_palm = surf_lsm_h%pavement_type(m)
3162          IF ( lup_palm == ind_lup_user )  THEN
3163             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3164             CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 )
3165          ELSEIF ( lup_palm == ind_lup_asph_conc )  THEN
3166             lup_dep = 9
3167          ELSEIF ( lup_palm == ind_lup_asph )  THEN
3168             lup_dep = 9
3169          ELSEIF ( lup_palm ==  ind_lup_conc )  THEN
3170             lup_dep = 9
3171          ELSEIF ( lup_palm ==  ind_lup_sett )  THEN
3172             lup_dep = 9
3173          ELSEIF ( lup_palm == ind_lup_pav_stones )  THEN
3174             lup_dep = 9
3175          ELSEIF ( lup_palm == ind_lup_cobblest )  THEN
3176             lup_dep = 9       
3177          ELSEIF ( lup_palm == ind_lup_metal )  THEN
3178             lup_dep = 9
3179          ELSEIF ( lup_palm == ind_lup_wood )  THEN
3180             lup_dep = 9   
3181          ELSEIF ( lup_palm == ind_lup_gravel )  THEN
3182             lup_dep = 9
3183          ELSEIF ( lup_palm == ind_lup_f_gravel )  THEN
3184             lup_dep = 9
3185          ELSEIF ( lup_palm == ind_lup_pebblest )  THEN
3186             lup_dep = 9
3187          ELSEIF ( lup_palm == ind_lup_woodchips )  THEN
3188             lup_dep = 9
3189          ELSEIF ( lup_palm == ind_lup_tartan )  THEN
3190             lup_dep = 9
3191          ELSEIF ( lup_palm == ind_lup_art_turf )  THEN
3192             lup_dep = 9
3193          ELSEIF ( lup_palm == ind_lup_clay )  THEN
3194             lup_dep = 9
3195          ENDIF
3196       ENDIF
3197
3198       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
3199          luw_palm = surf_lsm_h%water_type(m)     
3200          IF ( luw_palm == ind_luw_user )  THEN
3201             message_string = 'No water type defined. Please define water type to enable deposition calculation'
3202             CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 )
3203          ELSEIF ( luw_palm ==  ind_luw_lake )  THEN
3204             luw_dep = 13
3205          ELSEIF ( luw_palm == ind_luw_river )  THEN
3206             luw_dep = 13
3207          ELSEIF ( luw_palm == ind_luw_ocean )  THEN
3208             luw_dep = 6
3209          ELSEIF ( luw_palm == ind_luw_pond )  THEN
3210             luw_dep = 13
3211          ELSEIF ( luw_palm == ind_luw_fountain )  THEN
3212             luw_dep = 13 
3213          ENDIF
3214       ENDIF
3215
3216
3217       !
3218       !-- Set wetness indicator to dry or wet for lsm vegetation or pavement
3219       IF ( surf_lsm_h%c_liq(m) > 0 )  THEN
3220          nwet = 1
3221       ELSE
3222          nwet = 0
3223       ENDIF
3224
3225       !
3226       !--Compute length of time step
3227       IF ( call_chem_at_all_substeps )  THEN
3228          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
3229       ELSE
3230          dt_chem = dt_3d
3231       ENDIF
3232
3233
3234       dh = dzw(k)
3235       inv_dh = 1.0_wp / dh
3236       dt_dh = dt_chem/dh
3237
3238       !
3239       !-- Concentration at i,j,k
3240       DO  lsp = 1, nspec
3241          cc(lsp) = chem_species(lsp)%conc(k,j,i)
3242       ENDDO
3243
3244
3245       !
3246       !-- Temperature at i,j,k
3247       ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
3248
3249       ts       = ttemp - 273.15  !< in degrees celcius
3250
3251       !
3252       !-- Viscosity of air
3253       visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
3254
3255       !
3256       !-- Air density at k
3257       dens = rho_air_zw(k)
3258
3259       !
3260       !-- Calculate relative humidity from specific humidity for DEPAC
3261       qv_tmp = MAX(q(k,j,i),0.0_wp)
3262       relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) )
3263
3264
3265
3266       !
3267       !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
3268       !-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
3269
3270       !
3271       !-- Vegetation
3272       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
3273
3274
3275          slinnfac = 1.0_wp
3276
3277          !
3278          !-- Get vd
3279          DO  lsp = 1, nvar
3280             !
3281             !-- Initialize
3282             vs = 0.0_wp
3283             vd_lu = 0.0_wp
3284             Rs = 0.0_wp
3285             Rb = 0.0_wp
3286             Rc_tot = 0.0_wp
3287             IF ( spc_names(lsp) == 'PM10' )  THEN
3288                part_type = 1
3289                !
3290                !-- Sedimentation velocity
3291                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3292                     particle_pars(ind_p_size, part_type), &
3293                     particle_pars(ind_p_slip, part_type), &
3294                     visc)
3295
3296                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
3297                     vs, &
3298                     particle_pars(ind_p_size, part_type), &
3299                     particle_pars(ind_p_slip, part_type), &
3300                     nwet, ttemp, dens, visc, &
3301                     lu_dep,  &
3302                     r_aero_lu, ustar_lu )
3303
3304                bud_now_lu(lsp) = - cc(lsp) * &
3305                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3306
3307
3308             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3309                part_type = 2
3310                !
3311                !-- Sedimentation velocity
3312                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3313                     particle_pars(ind_p_size, part_type), &
3314                     particle_pars(ind_p_slip, part_type), &
3315                     visc)
3316
3317
3318                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
3319                     vs, &
3320                     particle_pars(ind_p_size, part_type), &
3321                     particle_pars(ind_p_slip, part_type), &
3322                     nwet, ttemp, dens, visc, &
3323                     lu_dep , &
3324                     r_aero_lu, ustar_lu )
3325
3326                bud_now_lu(lsp) = - cc(lsp) * &
3327                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3328
3329
3330             ELSE !< GASES
3331
3332                !
3333                !-- Read spc_name of current species for gas parameter
3334
3335                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
3336                   i_pspec = 0
3337                   DO  pspec = 1, nposp
3338                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3339                         i_pspec = pspec
3340                      END IF
3341                   ENDDO
3342
3343                ELSE
3344                   !
3345                   !-- For now species not deposited
3346                   CYCLE
3347                ENDIF
3348
3349                !
3350                !-- Factor used for conversion from ppb to ug/m3 :
3351                !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3352                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3353                !   c           1e-9              xm_tracer         1e9       /       xm_air            dens
3354                !-- thus:
3355                !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3356                !-- Use density at k:
3357
3358                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3359
3360                !
3361                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
3362                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3363                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
3364
3365                !
3366                !-- Diffusivity for DEPAC relevant gases
3367                !-- Use default value
3368                diffc            = 0.11e-4
3369                !
3370                !-- overwrite with known coefficients of diffusivity from Massman (1998)
3371                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
3372                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
3373                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
3374                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
3375                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
3376                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
3377                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
3378
3379
3380                !
3381                !-- Get quasi-laminar boundary layer resistance Rb:
3382                CALL get_rb_cell( (lu_dep == ilu_water_sea) .OR. (lu_dep == ilu_water_inland), &
3383                     z0h_lu, ustar_lu, diffc, &
3384                     Rb )
3385
3386                !
3387                !-- Get Rc_tot
3388                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
3389                     relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
3390                     r_aero_lu , Rb )
3391
3392
3393                !
3394                !-- Calculate budget
3395                IF ( Rc_tot <= 0.0 )  THEN
3396
3397                   bud_now_lu(lsp) = 0.0_wp
3398
3399                ELSE
3400
3401                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
3402
3403                   bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
3404                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3405                ENDIF
3406
3407             ENDIF
3408          ENDDO
3409       ENDIF
3410
3411
3412       !
3413       !-- Pavement
3414       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
3415
3416
3417          !
3418          !-- No vegetation on pavements:
3419          lai = 0.0_wp
3420          sai = 0.0_wp
3421
3422          slinnfac = 1.0_wp
3423
3424          !
3425          !-- Get vd
3426          DO  lsp = 1, nvar
3427             !
3428             !-- Initialize
3429             vs = 0.0_wp
3430             vd_lu = 0.0_wp
3431             Rs = 0.0_wp
3432             Rb = 0.0_wp
3433             Rc_tot = 0.0_wp
3434             IF ( spc_names(lsp) == 'PM10' )  THEN
3435                part_type = 1
3436                !
3437                !-- Sedimentation velocity:
3438                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3439                     particle_pars(ind_p_size, part_type), &
3440                     particle_pars(ind_p_slip, part_type), &
3441                     visc)
3442
3443                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
3444                     vs, &
3445                     particle_pars(ind_p_size, part_type), &
3446                     particle_pars(ind_p_slip, part_type), &
3447                     nwet, ttemp, dens, visc, &
3448                     lup_dep,  &
3449                     r_aero_lu, ustar_lu )
3450
3451                bud_now_lup(lsp) = - cc(lsp) * &
3452                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3453
3454
3455             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3456                part_type = 2
3457                !
3458                !-- Sedimentation velocity:
3459                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3460                     particle_pars(ind_p_size, part_type), &
3461                     particle_pars(ind_p_slip, part_type), &
3462                     visc)
3463
3464
3465                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
3466                     vs, &
3467                     particle_pars(ind_p_size, part_type), &
3468                     particle_pars(ind_p_slip, part_type), &
3469                     nwet, ttemp, dens, visc, &
3470                     lup_dep, &
3471                     r_aero_lu, ustar_lu )
3472
3473                bud_now_lup(lsp) = - cc(lsp) * &
3474                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3475
3476
3477             ELSE  !<GASES
3478
3479                !
3480                !-- Read spc_name of current species for gas parameter
3481
3482                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
3483                   i_pspec = 0
3484                   DO  pspec = 1, nposp
3485                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3486                         i_pspec = pspec
3487                      END IF
3488                   ENDDO
3489
3490                ELSE
3491                   !
3492                   !-- For now species not deposited
3493                   CYCLE
3494                ENDIF
3495
3496                !-- Factor used for conversion from ppb to ug/m3 :
3497                !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3498                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3499                !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
3500                !-- thus:
3501                !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3502                !-- Use density at lowest layer:
3503
3504                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3505
3506                !
3507                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
3508                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3509                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
3510
3511                !
3512                !-- Diffusivity for DEPAC relevant gases
3513                !-- Use default value
3514                diffc            = 0.11e-4
3515                !
3516                !-- overwrite with known coefficients of diffusivity from Massman (1998)
3517                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
3518                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
3519                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
3520                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
3521                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
3522                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
3523                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
3524
3525
3526                !
3527                !-- Get quasi-laminar boundary layer resistance Rb:
3528                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland), &
3529                     z0h_lu, ustar_lu, diffc, &
3530                     Rb )
3531
3532
3533                !
3534                !-- Get Rc_tot
3535                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
3536                     relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
3537                     r_aero_lu , Rb )
3538
3539
3540                !
3541                !-- Calculate budget
3542                IF ( Rc_tot <= 0.0 )  THEN
3543
3544                   bud_now_lup(lsp) = 0.0_wp
3545
3546                ELSE
3547
3548                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
3549
3550                   bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
3551                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3552                ENDIF
3553
3554
3555             ENDIF
3556          ENDDO
3557       ENDIF
3558
3559
3560       !
3561       !-- Water
3562       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
3563
3564
3565          !
3566          !-- No vegetation on water:
3567          lai = 0.0_wp
3568          sai = 0.0_wp
3569
3570          slinnfac = 1.0_wp
3571
3572          !
3573          !-- Get vd
3574          DO  lsp = 1, nvar
3575             !
3576             !-- Initialize
3577             vs = 0.0_wp
3578             vd_lu = 0.0_wp
3579             Rs = 0.0_wp
3580             Rb = 0.0_wp
3581             Rc_tot = 0.0_wp 
3582             IF ( spc_names(lsp) == 'PM10' )  THEN
3583                part_type = 1
3584                !
3585                !-- Sedimentation velocity:
3586                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3587                     particle_pars(ind_p_size, part_type), &
3588                     particle_pars(ind_p_slip, part_type), &
3589                     visc)
3590
3591                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
3592                     vs, &
3593                     particle_pars(ind_p_size, part_type), &
3594                     particle_pars(ind_p_slip, part_type), &
3595                     nwet, ttemp, dens, visc, &
3596                     luw_dep, &
3597                     r_aero_lu, ustar_lu )
3598
3599                bud_now_luw(lsp) = - cc(lsp) * &
3600                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3601
3602
3603             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3604                part_type = 2
3605                !
3606                !-- Sedimentation velocity:
3607                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3608                     particle_pars(ind_p_size, part_type), &
3609                     particle_pars(ind_p_slip, part_type), &
3610                     visc)
3611
3612
3613                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
3614                     vs, &
3615                     particle_pars(ind_p_size, part_type), &
3616                     particle_pars(ind_p_slip, part_type), &
3617                     nwet, ttemp, dens, visc, &
3618                     luw_dep, &
3619                     r_aero_lu, ustar_lu )
3620
3621                bud_now_luw(lsp) = - cc(lsp) * &
3622                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3623
3624
3625             ELSE  !<GASES
3626
3627                !
3628                !-- Read spc_name of current species for gas PARAMETER
3629
3630                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
3631                   i_pspec = 0
3632                   DO  pspec = 1, nposp
3633                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3634                         i_pspec = pspec
3635                      END IF
3636                   ENDDO
3637
3638                ELSE
3639                   !
3640                   !-- For now species not deposited
3641                   CYCLE
3642                ENDIF
3643
3644                !-- Factor used for conversion from ppb to ug/m3 :
3645                !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3646                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3647                !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
3648                !-- thus:
3649                !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3650                !-- Use density at lowest layer:
3651
3652                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3653
3654                !
3655                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
3656                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3657                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
3658
3659                !
3660                !-- Diffusivity for DEPAC relevant gases
3661                !-- Use default value
3662                diffc            = 0.11e-4
3663                !
3664                !-- overwrite with known coefficients of diffusivity from Massman (1998)
3665                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
3666                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
3667                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
3668                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
3669                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
3670                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
3671                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
3672
3673
3674                !
3675                !-- Get quasi-laminar boundary layer resistance Rb:
3676                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland), &
3677                     z0h_lu, ustar_lu, diffc, &
3678                     Rb )
3679
3680                !
3681                !-- Get Rc_tot
3682                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
3683                     relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
3684                     r_aero_lu , Rb )
3685
3686
3687                ! Calculate budget
3688                IF ( Rc_tot <= 0.0 )  THEN
3689
3690                   bud_now_luw(lsp) = 0.0_wp
3691
3692                ELSE
3693
3694                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
3695
3696                   bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
3697                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3698                ENDIF
3699
3700             ENDIF
3701          ENDDO
3702       ENDIF
3703
3704
3705       bud_now = 0.0_wp
3706       !
3707       !-- Calculate overall budget for surface m and adapt concentration
3708       DO  lsp = 1, nspec
3709
3710
3711          bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_now_lu(lsp) + &
3712               surf_lsm_h%frac(ind_pav_green,m) * bud_now_lup(lsp) + &
3713               surf_lsm_h%frac(ind_wat_win,m) * bud_now_luw(lsp)
3714
3715          !
3716          !-- Compute new concentration:
3717          cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
3718
3719          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, cc(lsp))
3720
3721       ENDDO
3722
3723    ENDIF
3724
3725
3726
3727
3728    !
3729    !-- For USM surfaces   
3730
3731    IF ( match_usm )  THEN
3732
3733       !
3734       !-- Get surface element information at i,j:
3735       m = surf_usm_h%start_index(j,i)
3736       k = surf_usm_h%k(m)
3737
3738       !
3739       !-- Get needed variables for surface element m
3740       ustar_lu  = surf_usm_h%us(m)
3741       z0h_lu    = surf_usm_h%z0h(m)
3742       r_aero_lu = surf_usm_h%r_a(m)
3743       glrad     = surf_usm_h%rad_sw_in(m)
3744       lai = surf_usm_h%lai(m)
3745       sai = lai + 1
3746
3747       !
3748       !-- For small grid spacing neglect R_a
3749       IF ( dzw(k) <= 1.0 )  THEN
3750          r_aero_lu = 0.0_wp
3751       ENDIF
3752
3753       !
3754       !-- Initialize lu's
3755       luu_palm = 0
3756       luu_dep = 0
3757       lug_palm = 0
3758       lug_dep = 0
3759       lud_palm = 0
3760       lud_dep = 0
3761
3762       !
3763       !-- Initialize budgets
3764       bud_now_luu  = 0.0_wp
3765       bud_now_lug = 0.0_wp
3766       bud_now_lud = 0.0_wp
3767
3768
3769       !
3770       !-- Get land use for i,j and assign to DEPAC lu
3771       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
3772          !
3773          !-- For green urban surfaces (e.g. green roofs
3774          !-- assume LU short grass
3775          lug_palm = ind_lu_s_grass
3776          IF ( lug_palm == ind_lu_user )  THEN
3777             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
3778             CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
3779          ELSEIF ( lug_palm == ind_lu_b_soil )  THEN
3780             lug_dep = 9
3781          ELSEIF ( lug_palm == ind_lu_mixed_crops )  THEN
3782             lug_dep = 2
3783          ELSEIF ( lug_palm == ind_lu_s_grass )  THEN
3784             lug_dep = 1
3785          ELSEIF ( lug_palm == ind_lu_ev_needle_trees )  THEN
3786             lug_dep = 4
3787          ELSEIF ( lug_palm == ind_lu_de_needle_trees )  THEN
3788             lug_dep = 4
3789          ELSEIF ( lug_palm == ind_lu_ev_broad_trees )  THEN
3790             lug_dep = 12
3791          ELSEIF ( lug_palm == ind_lu_de_broad_trees )  THEN
3792             lug_dep = 5
3793          ELSEIF ( lug_palm == ind_lu_t_grass )  THEN
3794             lug_dep = 1
3795          ELSEIF ( lug_palm == ind_lu_desert )  THEN
3796             lug_dep = 9
3797          ELSEIF ( lug_palm == ind_lu_tundra  )  THEN
3798             lug_dep = 8
3799          ELSEIF ( lug_palm == ind_lu_irr_crops )  THEN
3800             lug_dep = 2
3801          ELSEIF ( lug_palm == ind_lu_semidesert )  THEN
3802             lug_dep = 8
3803          ELSEIF ( lug_palm == ind_lu_ice )  THEN
3804             lug_dep = 10
3805          ELSEIF ( lug_palm == ind_lu_marsh )  THEN
3806             lug_dep = 8
3807          ELSEIF ( lug_palm == ind_lu_ev_shrubs )  THEN
3808             lug_dep = 14
3809          ELSEIF ( lug_palm == ind_lu_de_shrubs  )  THEN
3810             lug_dep = 14
3811          ELSEIF ( lug_palm == ind_lu_mixed_forest )  THEN
3812             lug_dep = 4
3813          ELSEIF ( lug_palm == ind_lu_intrup_forest )  THEN
3814             lug_dep = 8     
3815          ENDIF
3816       ENDIF
3817
3818       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
3819          !
3820          !-- For walls in USM assume concrete walls/roofs,
3821          !-- assumed LU class desert as also assumed for
3822          !-- pavements in LSM
3823          luu_palm = ind_lup_conc
3824          IF ( luu_palm == ind_lup_user )  THEN
3825             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3826             CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 )
3827          ELSEIF ( luu_palm == ind_lup_asph_conc )  THEN
3828             luu_dep = 9
3829          ELSEIF ( luu_palm == ind_lup_asph )  THEN
3830             luu_dep = 9
3831          ELSEIF ( luu_palm ==  ind_lup_conc )  THEN
3832             luu_dep = 9
3833          ELSEIF ( luu_palm ==  ind_lup_sett )  THEN
3834             luu_dep = 9
3835          ELSEIF ( luu_palm == ind_lup_pav_stones )  THEN
3836             luu_dep = 9
3837          ELSEIF ( luu_palm == ind_lup_cobblest )  THEN
3838             luu_dep = 9       
3839          ELSEIF ( luu_palm == ind_lup_metal )  THEN
3840             luu_dep = 9
3841          ELSEIF ( luu_palm == ind_lup_wood )  THEN
3842             luu_dep = 9   
3843          ELSEIF ( luu_palm == ind_lup_gravel )  THEN
3844             luu_dep = 9
3845          ELSEIF ( luu_palm == ind_lup_f_gravel )  THEN
3846             luu_dep = 9
3847          ELSEIF ( luu_palm == ind_lup_pebblest )  THEN
3848             luu_dep = 9
3849          ELSEIF ( luu_palm == ind_lup_woodchips )  THEN
3850             luu_dep = 9
3851          ELSEIF ( luu_palm == ind_lup_tartan )  THEN
3852             luu_dep = 9
3853          ELSEIF ( luu_palm == ind_lup_art_turf )  THEN
3854             luu_dep = 9
3855          ELSEIF ( luu_palm == ind_lup_clay )  THEN
3856             luu_dep = 9
3857          ENDIF
3858       ENDIF
3859
3860       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
3861          !
3862          !-- For windows in USM assume metal as this is
3863          !-- as close as we get, assumed LU class desert
3864          !-- as also assumed for pavements in LSM
3865          lud_palm = ind_lup_metal     
3866          IF ( lud_palm == ind_lup_user )  THEN
3867             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3868             CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 )
3869          ELSEIF ( lud_palm == ind_lup_asph_conc )  THEN
3870             lud_dep = 9
3871          ELSEIF ( lud_palm == ind_lup_asph )  THEN
3872             lud_dep = 9
3873          ELSEIF ( lud_palm ==  ind_lup_conc )  THEN
3874             lud_dep = 9
3875          ELSEIF ( lud_palm ==  ind_lup_sett )  THEN
3876             lud_dep = 9
3877          ELSEIF ( lud_palm == ind_lup_pav_stones )  THEN
3878             lud_dep = 9
3879          ELSEIF ( lud_palm == ind_lup_cobblest )  THEN
3880             lud_dep = 9       
3881          ELSEIF ( lud_palm == ind_lup_metal )  THEN
3882             lud_dep = 9
3883          ELSEIF ( lud_palm == ind_lup_wood )  THEN
3884             lud_dep = 9   
3885          ELSEIF ( lud_palm == ind_lup_gravel )  THEN
3886             lud_dep = 9
3887          ELSEIF ( lud_palm == ind_lup_f_gravel )  THEN
3888             lud_dep = 9
3889          ELSEIF ( lud_palm == ind_lup_pebblest )  THEN
3890             lud_dep = 9
3891          ELSEIF ( lud_palm == ind_lup_woodchips )  THEN
3892             lud_dep = 9
3893          ELSEIF ( lud_palm == ind_lup_tartan )  THEN
3894             lud_dep = 9
3895          ELSEIF ( lud_palm == ind_lup_art_turf )  THEN
3896             lud_dep = 9
3897          ELSEIF ( lud_palm == ind_lup_clay )  THEN
3898             lud_dep = 9
3899          ENDIF
3900       ENDIF
3901
3902
3903       !
3904       !-- @TODO: Activate these lines as soon as new ebsolver branch is merged:
3905       !-- Set wetness indicator to dry or wet for usm vegetation or pavement
3906       !IF ( surf_usm_h%c_liq(m) > 0 )  THEN
3907       !   nwet = 1
3908       !ELSE
3909       nwet = 0
3910       !ENDIF
3911
3912       !
3913       !-- Compute length of time step
3914       IF ( call_chem_at_all_substeps )  THEN
3915          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
3916       ELSE
3917          dt_chem = dt_3d
3918       ENDIF
3919
3920
3921       dh = dzw(k)
3922       inv_dh = 1.0_wp / dh
3923       dt_dh = dt_chem/dh
3924
3925       !
3926       !-- Concentration at i,j,k
3927       DO  lsp = 1, nspec
3928          cc(lsp) = chem_species(lsp)%conc(k,j,i)
3929       ENDDO
3930
3931       !
3932       !-- Temperature at i,j,k
3933       ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
3934
3935       ts       = ttemp - 273.15  !< in degrees celcius
3936
3937       !
3938       !-- Viscosity of air
3939       visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
3940
3941       !
3942       !-- Air density at k
3943       dens = rho_air_zw(k)
3944
3945       !
3946       !-- Calculate relative humidity from specific humidity for DEPAC
3947       qv_tmp = MAX(q(k,j,i),0.0_wp)
3948       relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(nzb) )
3949
3950
3951
3952       !
3953       !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
3954       !-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
3955
3956       !
3957       !-- Walls/roofs
3958       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
3959
3960
3961          !
3962          !-- No vegetation on non-green walls:
3963          lai = 0.0_wp
3964          sai = 0.0_wp
3965
3966          slinnfac = 1.0_wp
3967
3968          !
3969          !-- Get vd
3970          DO  lsp = 1, nvar
3971             !
3972             !-- Initialize
3973             vs = 0.0_wp
3974             vd_lu = 0.0_wp
3975             Rs = 0.0_wp
3976             Rb = 0.0_wp
3977             Rc_tot = 0.0_wp
3978             IF (spc_names(lsp) == 'PM10' )  THEN
3979                part_type = 1
3980                !
3981                !-- Sedimentation velocity
3982                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3983                     particle_pars(ind_p_size, part_type), &
3984                     particle_pars(ind_p_slip, part_type), &
3985                     visc)
3986
3987                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
3988                     vs, &
3989                     particle_pars(ind_p_size, part_type), &
3990                     particle_pars(ind_p_slip, part_type), &
3991                     nwet, ttemp, dens, visc, &
3992                     luu_dep,  &
3993                     r_aero_lu, ustar_lu )
3994
3995                bud_now_luu(lsp) = - cc(lsp) * &
3996                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3997
3998
3999             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4000                part_type = 2
4001                !
4002                !-- Sedimentation velocity
4003                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4004                     particle_pars(ind_p_size, part_type), &
4005                     particle_pars(ind_p_slip, part_type), &
4006                     visc)
4007
4008
4009                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
4010                     vs, &
4011                     particle_pars(ind_p_size, part_type), &
4012                     particle_pars(ind_p_slip, part_type), &
4013                     nwet, ttemp, dens, visc, &
4014                     luu_dep , &
4015                     r_aero_lu, ustar_lu )
4016
4017                bud_now_luu(lsp) = - cc(lsp) * &
4018                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4019
4020
4021             ELSE  !< GASES
4022
4023                !
4024                !-- Read spc_name of current species for gas parameter
4025
4026                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4027                   i_pspec = 0
4028                   DO  pspec = 1, nposp
4029                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4030                         i_pspec = pspec
4031                      END IF
4032                   ENDDO
4033                ELSE
4034                   !
4035                   !-- For now species not deposited
4036                   CYCLE
4037                ENDIF
4038
4039                !
4040                !-- Factor used for conversion from ppb to ug/m3 :
4041                !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4042                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4043                !   c           1e-9              xm_tracer         1e9       /       xm_air            dens
4044                !-- thus:
4045                !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4046                !-- Use density at k:
4047
4048                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
4049
4050                !
4051                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
4052                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4053                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
4054
4055                !
4056                !-- Diffusivity for DEPAC relevant gases
4057                !-- Use default value
4058                diffc            = 0.11e-4
4059                !
4060                !-- overwrite with known coefficients of diffusivity from Massman (1998)
4061                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
4062                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
4063                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
4064                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
4065                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
4066                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
4067                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
4068
4069
4070                !
4071                !-- Get quasi-laminar boundary layer resistance Rb:
4072                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland), &
4073                     z0h_lu, ustar_lu, diffc, &
4074                     Rb )
4075
4076                !
4077                !-- Get Rc_tot
4078                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
4079                     relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
4080                     r_aero_lu , Rb )
4081
4082
4083                !
4084                !-- Calculate budget
4085                IF ( Rc_tot <= 0.0 )  THEN
4086
4087                   bud_now_luu(lsp) = 0.0_wp
4088
4089                ELSE
4090
4091                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
4092
4093                   bud_now_luu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
4094                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4095                ENDIF
4096
4097             ENDIF
4098          ENDDO
4099       ENDIF
4100
4101
4102       !
4103       !-- Green usm surfaces
4104       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
4105
4106
4107          slinnfac = 1.0_wp
4108
4109          !
4110          !-- Get vd
4111          DO  lsp = 1, nvar
4112             !
4113             !-- Initialize
4114             vs = 0.0_wp
4115             vd_lu = 0.0_wp
4116             Rs = 0.0_wp
4117             Rb = 0.0_wp
4118             Rc_tot = 0.0_wp
4119             IF ( spc_names(lsp) == 'PM10' )  THEN
4120                part_type = 1
4121                ! Sedimentation velocity
4122                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4123                     particle_pars(ind_p_size, part_type), &
4124                     particle_pars(ind_p_slip, part_type), &
4125                     visc)
4126
4127                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
4128                     vs, &
4129                     particle_pars(ind_p_size, part_type), &
4130                     particle_pars(ind_p_slip, part_type), &
4131                     nwet, ttemp, dens, visc, &
4132                     lug_dep,  &
4133                     r_aero_lu, ustar_lu )
4134
4135                bud_now_lug(lsp) = - cc(lsp) * &
4136                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4137
4138
4139             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4140                part_type = 2
4141                ! Sedimentation velocity
4142                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4143                     particle_pars(ind_p_size, part_type), &
4144                     particle_pars(ind_p_slip, part_type), &
4145                     visc)
4146
4147
4148                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
4149                     vs, &
4150                     particle_pars(ind_p_size, part_type), &
4151                     particle_pars(ind_p_slip, part_type), &
4152                     nwet, ttemp, dens, visc, &
4153                     lug_dep, &
4154                     r_aero_lu, ustar_lu )
4155
4156                bud_now_lug(lsp) = - cc(lsp) * &
4157                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4158
4159
4160             ELSE  !< GASES
4161
4162                !
4163                !-- Read spc_name of current species for gas parameter
4164
4165                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4166                   i_pspec = 0
4167                   DO  pspec = 1, nposp
4168                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4169                         i_pspec = pspec
4170                      END IF
4171                   ENDDO
4172                ELSE
4173                   !
4174                   !-- For now species not deposited
4175                   CYCLE
4176                ENDIF
4177
4178                !
4179                !-- Factor used for conversion from ppb to ug/m3 :
4180                !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4181                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4182                !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
4183                !-- thus:
4184                !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4185                !-- Use density at k:
4186
4187                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
4188
4189                !
4190                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
4191                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4192                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
4193
4194                !
4195                !-- Diffusivity for DEPAC relevant gases
4196                !-- Use default value
4197                diffc            = 0.11e-4
4198                !
4199                !-- overwrite with known coefficients of diffusivity from Massman (1998)
4200                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
4201                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
4202                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
4203                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
4204                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
4205                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
4206                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
4207
4208
4209                !
4210                !-- Get quasi-laminar boundary layer resistance Rb:
4211                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland), &
4212                     z0h_lu, ustar_lu, diffc, &
4213                     Rb )
4214
4215
4216                !
4217                !-- Get Rc_tot
4218                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
4219                     relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
4220                     r_aero_lu , Rb )
4221
4222
4223                !
4224                !-- Calculate budget
4225                IF ( Rc_tot <= 0.0 )  THEN
4226
4227                   bud_now_lug(lsp) = 0.0_wp
4228
4229                ELSE
4230
4231                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
4232
4233                   bud_now_lug(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
4234                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4235                ENDIF
4236
4237
4238             ENDIF
4239          ENDDO
4240       ENDIF
4241
4242
4243       !
4244       !-- Windows
4245       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
4246
4247
4248          !
4249          !-- No vegetation on windows:
4250          lai = 0.0_wp
4251          sai = 0.0_wp
4252
4253          slinnfac = 1.0_wp
4254
4255          !
4256          !-- Get vd
4257          DO  lsp = 1, nvar
4258             !
4259             !-- Initialize
4260             vs = 0.0_wp
4261             vd_lu = 0.0_wp
4262             Rs = 0.0_wp
4263             Rb = 0.0_wp
4264             Rc_tot = 0.0_wp 
4265             IF ( spc_names(lsp) == 'PM10' )  THEN
4266                part_type = 1
4267                !
4268                !-- Sedimentation velocity
4269                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4270                     particle_pars(ind_p_size, part_type), &
4271                     particle_pars(ind_p_slip, part_type), &
4272                     visc)
4273
4274                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
4275                     vs, &
4276                     particle_pars(ind_p_size, part_type), &
4277                     particle_pars(ind_p_slip, part_type), &
4278                     nwet, ttemp, dens, visc, &
4279                     lud_dep, &
4280                     r_aero_lu, ustar_lu )
4281
4282                bud_now_lud(lsp) = - cc(lsp) * &
4283                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4284
4285
4286             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4287                part_type = 2
4288                !
4289                !-- Sedimentation velocity
4290                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4291                     particle_pars(ind_p_size, part_type), &
4292                     particle_pars(ind_p_slip, part_type), &
4293                     visc)
4294
4295
4296                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
4297                     vs, &
4298                     particle_pars(ind_p_size, part_type), &
4299                     particle_pars(ind_p_slip, part_type), &
4300                     nwet, ttemp, dens, visc, &
4301                     lud_dep, &
4302                     r_aero_lu, ustar_lu )
4303
4304                bud_now_lud(lsp) = - cc(lsp) * &
4305                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4306
4307
4308             ELSE  !< GASES
4309
4310                !
4311                !-- Read spc_name of current species for gas PARAMETER
4312
4313                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4314                   i_pspec = 0
4315                   DO  pspec = 1, nposp
4316                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4317                         i_pspec = pspec
4318                      END IF
4319                   ENDDO
4320                ELSE
4321                   ! For now species not deposited
4322                   CYCLE
4323                ENDIF
4324
4325                !
4326                !-- Factor used for conversion from ppb to ug/m3 :
4327                !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4328                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4329                !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
4330                !-- thus:
4331                !    c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4332                !-- Use density at k:
4333
4334                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
4335
4336                !
4337                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
4338                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4339                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
4340
4341                !
4342                !-- Diffusivity for DEPAC relevant gases
4343                !-- Use default value
4344                diffc            = 0.11e-4
4345                !
4346                !-- overwrite with known coefficients of diffusivity from Massman (1998)
4347                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4 
4348                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
4349                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
4350                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
4351                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
4352                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
4353                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
4354
4355
4356                !
4357                !-- Get quasi-laminar boundary layer resistance Rb:
4358                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland), &
4359                     z0h_lu, ustar_lu, diffc, &
4360                     Rb )
4361
4362                !
4363                !-- Get Rc_tot
4364                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
4365                     relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
4366                     r_aero_lu , Rb )
4367
4368
4369                !
4370                !-- Calculate budget
4371                IF ( Rc_tot <= 0.0 )  THEN
4372
4373                   bud_now_lud(lsp) = 0.0_wp
4374
4375                ELSE
4376
4377                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
4378
4379                   bud_now_lud(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
4380                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4381                ENDIF
4382
4383             ENDIF
4384          ENDDO
4385       ENDIF
4386
4387
4388       bud_now = 0.0_wp
4389       !
4390       !-- Calculate overall budget for surface m and adapt concentration
4391       DO  lsp = 1, nspec
4392
4393
4394          bud_now(lsp) = surf_usm_h%frac(ind_veg_wall,m) * bud_now_luu(lsp) + &
4395               surf_usm_h%frac(ind_pav_green,m) * bud_now_lug(lsp) + &
4396               surf_usm_h%frac(ind_wat_win,m) * bud_now_lud(lsp)
4397
4398          !
4399          !-- Compute new concentration
4400          cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
4401
4402          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, cc(lsp) )
4403
4404       ENDDO
4405
4406    ENDIF
4407
4408
4409
4410 END SUBROUTINE chem_depo
4411
4412
4413
4414 !----------------------------------------------------------------------------------
4415 !
4416 !> DEPAC:
4417 !> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
4418 !> module of the LOTOS-EUROS model (Manders et al., 2017)
4419 !   
4420 !> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
4421 !> van Zanten et al., 2010.
4422 !---------------------------------------------------------------------------------
4423 !   
4424 !----------------------------------------------------------------------------------   
4425 !
4426 !> depac: compute total canopy (or surface) resistance Rc for gases
4427 !----------------------------------------------------------------------------------
4428
4429 SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, glrad, sinphi,  &
4430      rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, catm, diffc,           &
4431      ra, rb ) 
4432
4433
4434   !
4435   !--Some of depac arguments are OPTIONAL:
4436   !
4437   !-- A. compute Rc_tot without compensation points (ccomp_tot will be zero):
4438   !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
4439   !
4440   !-- B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
4441   !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4442   !                  c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
4443   !
4444   !-- C. compute effective Rc based on compensation points (used for OPS):
4445   !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4446   !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4447   !                 ra, rb, rc_eff)
4448   !-- X1. Extra (OPTIONAL) output variables:
4449   !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4450   !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4451   !                 ra, rb, rc_eff, &
4452   !                 gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
4453   !-- X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
4454   !--     CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4455   !                 c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4456   !                 ra, rb, rc_eff, &
4457   !                 gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
4458   !                 calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
4459
4460
4461    IMPLICIT NONE
4462
4463    CHARACTER(len=*), INTENT(IN) ::  compnam         !< component name
4464                                                     !< 'HNO3','NO','NO2','O3','SO2','NH3'
4465    INTEGER(iwp), INTENT(IN) ::  day_of_year         !< day of year, 1 ... 365 (366)
4466    INTEGER(iwp), INTENT(IN) ::  nwet                !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4467    INTEGER(iwp), INTENT(IN) ::  lu                  !< land use type, lu = 1,...,nlu
4468    INTEGER(iwp), INTENT(IN) ::  iratns              !< index for NH3/SO2 ratio used for SO2:
4469                                                     !< iratns = 1: low NH3/SO2
4470                                                     !< iratns = 2: high NH3/SO2
4471                                                     !< iratns = 3: very low NH3/SO2
4472    REAL(wp), INTENT(IN) ::  lat                     !< latitude Northern hemisphere (degrees) (S. hemisphere not possible)
4473    REAL(wp), INTENT(IN) ::  t                       !< temperature (C)
4474    REAL(wp), INTENT(IN) ::  ust                     !< friction velocity (m/s)
4475    REAL(wp), INTENT(IN) ::  glrad                   !< global radiation (W/m2)
4476    REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
4477    REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
4478    REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
4479    REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and stems)
4480    REAL(wp), INTENT(IN) ::  diffc                   !< diffusivity
4481    REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
4482    REAL(wp), INTENT(IN) ::  catm                    !< actual atmospheric concentration (ug/m3)
4483    REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
4484    REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
4485
4486    REAL(wp), INTENT(OUT) ::  rc_tot                 !< total canopy resistance Rc (s/m)
4487    REAL(wp), INTENT(OUT) ::  ccomp_tot              !< total compensation point (ug/m3)
4488                                                     !< [= 0 for species that don't have a compensation
4489
4490    !
4491    !-- Local variables:
4492    !
4493    !-- Component number taken from component name, paramteres matched with include files
4494    INTEGER(iwp) ::  icmp
4495
4496    !
4497    !-- Component numbers:
4498    INTEGER(iwp), PARAMETER ::  icmp_o3   = 1
4499    INTEGER(iwp), PARAMETER ::  icmp_so2  = 2
4500    INTEGER(iwp), PARAMETER ::  icmp_no2  = 3
4501    INTEGER(iwp), PARAMETER ::  icmp_no   = 4
4502    INTEGER(iwp), PARAMETER ::  icmp_nh3  = 5
4503    INTEGER(iwp), PARAMETER ::  icmp_co   = 6
4504    INTEGER(iwp), PARAMETER ::  icmp_no3  = 7
4505    INTEGER(iwp), PARAMETER ::  icmp_hno3 = 8
4506    INTEGER(iwp), PARAMETER ::  icmp_n2o5 = 9
4507    INTEGER(iwp), PARAMETER ::  icmp_h2o2 = 10
4508
4509    LOGICAL ::  ready                                !< Rc has been set:
4510    !< = 1 -> constant Rc
4511    !< = 2 -> temperature dependent Rc
4512    !
4513    !-- Vegetation indicators:
4514    LOGICAL ::  lai_present                          !< leaves are present for current land use type
4515    LOGICAL ::  sai_present                          !< vegetation is present for current land use type
4516
4517    REAL(wp) ::  laimax                              !< maximum leaf area index (-)
4518    REAL(wp) ::  gw                                  !< external leaf conductance (m/s)
4519    REAL(wp) ::  gstom                               !< stomatal conductance (m/s)
4520    REAL(wp) ::  gsoil_eff                           !< effective soil conductance (m/s)
4521    REAL(wp) ::  gc_tot                              !< total canopy conductance (m/s)
4522    REAL(wp) ::  cw                                  !< external leaf surface compensation point (ug/m3)
4523    REAL(wp) ::  cstom                               !< stomatal compensation point (ug/m3)
4524    REAL(wp) ::  csoil                               !< soil compensation point (ug/m3)
4525
4526
4527
4528    !
4529    !-- Define component number
4530    SELECT CASE ( TRIM( compnam ) )
4531
4532    CASE ( 'O3', 'o3' )
4533       icmp = icmp_o3
4534
4535    CASE ( 'SO2', 'so2' )
4536       icmp = icmp_so2
4537
4538    CASE ( 'NO2', 'no2' )
4539       icmp = icmp_no2
4540
4541    CASE ( 'NO', 'no' )
4542       icmp = icmp_no 
4543
4544    CASE ( 'NH3', 'nh3' )
4545       icmp = icmp_nh3
4546
4547    CASE ( 'CO', 'co' )
4548       icmp = icmp_co
4549
4550    CASE ( 'NO3', 'no3' )
4551       icmp = icmp_no3
4552
4553    CASE ( 'HNO3', 'hno3' )
4554       icmp = icmp_hno3
4555
4556    CASE ( 'N2O5', 'n2o5' )
4557       icmp = icmp_n2o5
4558
4559    CASE ( 'H2O2', 'h2o2' )
4560       icmp = icmp_h2o2
4561
4562    CASE default
4563       !
4564       !-- Component not part of DEPAC --> not deposited
4565       RETURN
4566
4567    END SELECT
4568
4569    !
4570    !-- Inititalize
4571    gw        = 0.0_wp
4572    gstom     = 0.0_wp
4573    gsoil_eff = 0.0_wp
4574    gc_tot    = 0.0_wp
4575    cw        = 0.0_wp
4576    cstom     = 0.0_wp
4577    csoil     = 0.0_wp
4578
4579
4580    !
4581    !-- Check whether vegetation is present:
4582    lai_present = ( lai > 0.0 )
4583    sai_present = ( sai > 0.0 )
4584
4585    !
4586    !-- Set Rc (i.e. rc_tot) in special cases:
4587    CALL rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
4588
4589    !
4590    !-- If Rc is not set:
4591    IF ( .NOT. ready ) then
4592
4593       !
4594       !-- External conductance:
4595       CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )         
4596
4597       !
4598       !-- Stomatal conductance:
4599       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, glrad, sinphi, t, rh, diffc, gstom, p )
4600       !
4601       !-- Effective soil conductance:
4602       CALL rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff )
4603
4604       !
4605       !-- Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
4606       CALL rc_rctot( gstom, gsoil_eff, gw, gc_tot, rc_tot )
4607
4608       !
4609       !-- Needed to include compensation point for NH3
4610       !-- Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
4611       !-- CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
4612       !     lai_present, sai_present, &
4613       !     ccomp_tot, &
4614       !     catm=catm,c_ave_prev_nh3=c_ave_prev_nh3, &
4615       !     c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
4616       !     tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
4617       !
4618       !-- Effective Rc based on compensation points:
4619       !   IF ( present(rc_eff) ) then
4620       !     check on required arguments:
4621       !      IF ( (.not. present(catm)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
4622       !         stop 'output argument rc_eff requires input arguments catm, ra and rb'
4623       !      END IF
4624       !--    Compute rc_eff :
4625       !      CALL rc_comp_point_rc_eff(ccomp_tot,catm,ra,rb,rc_tot,rc_eff)
4626       !   ENDIF
4627    ENDIF
4628
4629 END SUBROUTINE drydepos_gas_depac
4630
4631
4632
4633 !-------------------------------------------------------------------
4634 !> rc_special: compute total canopy resistance in special CASEs
4635 !-------------------------------------------------------------------
4636 SUBROUTINE rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
4637
4638    IMPLICIT NONE
4639   
4640    CHARACTER(len=*), INTENT(IN) ::  compnam     !< component name
4641
4642    INTEGER(iwp), INTENT(IN) ::  icmp            !< component index     
4643    INTEGER(iwp), INTENT(IN) ::  lu              !< land use type, lu = 1,...,nlu
4644    INTEGER(iwp), INTENT(IN) ::  nwet            !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4645
4646    REAL(wp), INTENT(IN) ::  t                   !< temperature (C)
4647
4648    REAL(wp), INTENT(OUT) ::  rc_tot             !< total canopy resistance Rc (s/m)
4649    REAL(wp), INTENT(OUT) ::  ccomp_tot          !< total compensation point (ug/m3)
4650
4651    LOGICAL, INTENT(OUT) ::  ready               !< Rc has been set
4652                                                 !< = 1 -> constant Rc
4653                                                 !< = 2 -> temperature dependent Rc
4654
4655    !
4656    !-- rc_tot is not yet set:
4657    ready = .FALSE.
4658
4659    !
4660    !-- Default compensation point in special CASEs = 0:
4661    ccomp_tot = 0.0_wp
4662
4663    SELECT CASE( TRIM( compnam ) )
4664    CASE( 'HNO3', 'N2O5', 'NO3', 'H2O2' )
4665       !
4666       !-- No separate resistances for HNO3; just one total canopy resistance:
4667       IF ( t < -5.0_wp .AND. nwet == 9 )  THEN
4668          !
4669          !-- T < 5 C and snow:
4670          rc_tot = 50.0_wp
4671       ELSE
4672          !
4673          !-- all other circumstances:
4674          rc_tot = 10.0_wp
4675       ENDIF
4676       ready = .TRUE.
4677
4678    CASE( 'NO', 'CO' )
4679       IF ( lu == ilu_water_sea .OR. lu == ilu_water_inland )  THEN       ! water
4680          rc_tot = 2000.0_wp
4681          ready = .TRUE.
4682       ELSEIF ( nwet == 1 )  THEN  !< wet
4683          rc_tot = 2000.0_wp
4684          ready = .TRUE.
4685       ENDIF
4686    CASE( 'NO2', 'O3', 'SO2', 'NH3' )
4687       !
4688       !-- snow surface:
4689       IF ( nwet == 9 )  THEN
4690          !
4691          !-- To be activated when snow is implemented
4692          !CALL rc_snow(ipar_snow(icmp),t,rc_tot)
4693          ready = .TRUE.
4694       ENDIF
4695    CASE default
4696       message_string = 'Component '// TRIM( compnam ) // ' not supported'
4697       CALL message( 'rc_special', 'CM0457', 1, 2, 0, 6, 0 )
4698    END SELECT
4699
4700 END SUBROUTINE rc_special
4701
4702
4703
4704 !-------------------------------------------------------------------
4705 !> rc_gw: compute external conductance
4706 !-------------------------------------------------------------------
4707 SUBROUTINE rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai, gw )
4708
4709    IMPLICIT NONE
4710
4711    !
4712    !-- Input/output variables:
4713    CHARACTER(len=*), INTENT(IN) ::  compnam      !< component name
4714
4715    INTEGER(iwp), INTENT(IN) ::  nwet             !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4716    INTEGER(iwp), INTENT(IN) ::  iratns           !< index for NH3/SO2 ratio;
4717                                                  !< iratns = 1: low NH3/SO2
4718                                                  !< iratns = 2: high NH3/SO2
4719                                                  !< iratns = 3: very low NH3/SO2
4720
4721    LOGICAL, INTENT(IN) ::  sai_present
4722
4723    REAL(wp), INTENT(IN) ::  t                    !< temperature (C)
4724    REAL(wp), INTENT(IN) ::  rh                   !< relative humidity (%)
4725    REAL(wp<