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

Last change on this file since 3833 was 3833, checked in by forkel, 2 years ago

removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead

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