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

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

Changed character length of name in species_def and photols_def to 15

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