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

Last change on this file since 3637 was 3637, checked in by knoop, 6 years ago

M Makefile

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