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

Last change on this file since 3458 was 3458, checked in by kanani, 4 years ago

Reintegrated fixes/changes from branch chemistry

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