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

Last change on this file since 3524 was 3524, checked in by raasch, 4 years ago

unused variables removed, missing working precision added, missing preprocessor directives added, bugfix concerning allocation of t_surf_wall_v in nopointer case, declaration statements rearranged to avoid compile time errors, mpi_abort arguments replaced to avoid compile errors

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