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

Last change on this file since 3611 was 3611, checked in by banzhafs, 3 years ago

chem_emissions_mod and chem_modules update to comply PALM coding rules

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