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

Last change on this file since 4004 was 4004, checked in by suehring, 2 years ago

chemistry: perform basic checks only when anthropenic emissions are switched on; virtual flights: allow arbitrary start/end positions also in return mode; bugfix in 2d data output

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