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

Last change on this file since 3700 was 3700, checked in by knoop, 4 years ago

Moved user_define_netdf_grid into user_module.f90
Added module interface for the definition of additional timeseries

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