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

Last change on this file since 3570 was 3570, checked in by kanani, 5 years ago

Fix too long lines (chemistry_model_mod, chem_emissions_mod), correct terminal message (palmrun)

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