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

Last change on this file since 3737 was 3737, checked in by suehring, 3 years ago

Enable initialization of chemistry variables via dynamic input file; enable mesoscale offline nesting of chemistry variables if data is in dynamic input file

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