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

Last change on this file since 3638 was 3638, checked in by forkel, 5 years ago

chemistry_model_mod: Added missing conversion factor fr2ppm for qvap

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