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

Last change on this file since 4069 was 4069, checked in by Giersch, 2 years ago

Bugfix for masked output, compiler warning removed, test case for wind turbine model revised

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