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

Last change on this file since 3652 was 3652, checked in by forkel, 5 years ago

Checks added for chemistry mechanism, parameter chem_mechanism added

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