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

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

Bugfix: set found logical correct in chem_data_output_2d

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