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

Last change on this file since 3820 was 3820, checked in by forkel, 3 years ago

renaming of get_mechanismname, do_emiss and do_depo, sorting in chem_modules

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