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

Last change on this file since 4370 was 4370, checked in by raasch, 19 months ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

  • Property svn:keywords set to Id
File size: 236.5 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-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2019 Karlsruhe Institute of Technology
19! Copyright 2017-2019 Freie Universitaet Berlin
20!------------------------------------------------------------------------------!
21!
22! Current revisions:
23! -----------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: chemistry_model_mod.f90 4370 2020-01-10 14:00:44Z raasch $
29! vector directives added to force vectorization on Intel19 compiler
30!
31! 4346 2019-12-18 11:55:56Z motisi
32! Introduction of wall_flags_total_0, which currently sets bits based on static
33! topography information used in wall_flags_static_0
34!
35! 4329 2019-12-10 15:46:36Z motisi
36! Renamed wall_flags_0 to wall_flags_static_0
37!
38! 4306 2019-11-25 12:04:48Z banzhafs
39! Corretion for r4304 commit
40!
41! 4304 2019-11-25 10:43:03Z banzhafs
42! Precision clean-up in drydepo_aero_zhang_vd subroutine
43!
44! 4292 2019-11-11 13:04:50Z banzhafs
45! Bugfix for r4290
46!
47! 4290 2019-11-11 12:06:14Z banzhafs
48! Bugfix in sedimentation resistance calculation in drydepo_aero_zhang_vd subroutine
49!
50! 4273 2019-10-24 13:40:54Z monakurppa
51! Add logical switches nesting_chem and nesting_offline_chem (both .TRUE.
52! by default)
53!
54! 4272 2019-10-23 15:18:57Z schwenkel
55! Further modularization of boundary conditions: moved boundary conditions to
56! respective modules
57!
58! 4268 2019-10-17 11:29:38Z schwenkel
59! Moving module specific boundary conditions from time_integration to module
60!
61! 4230 2019-09-11 13:58:14Z suehring
62! Bugfix, initialize mean profiles also in restart runs. Also initialize
63! array used for Runge-Kutta tendecies in restart runs. 
64!
65! 4227 2019-09-10 18:04:34Z gronemeier
66! implement new palm_date_time_mod
67!
68! 4182 2019-08-22 15:20:23Z scharf
69! Corrected "Former revisions" section
70!
71! 4167 2019-08-16 11:01:48Z suehring
72! Changed behaviour of masked output over surface to follow terrain and ignore
73! buildings (J.Resler, T.Gronemeier)
74!
75!
76! 4166 2019-08-16 07:54:21Z resler
77! Bugfix in decycling
78!
79! 4115 2019-07-24 12:50:49Z suehring
80! Fix faulty IF statement in decycling initialization
81!
82! 4110 2019-07-22 17:05:21Z suehring
83! - Decycling boundary conditions are only set at the ghost points not on the
84!   prognostic grid points
85! - Allocation and initialization of special advection flags cs_advc_flags_s
86!   used for chemistry. These are exclusively used for chemical species in
87!   order to distinguish from the usually-used flags which might be different
88!   when decycling is applied in combination with cyclic boundary conditions.
89!   Moreover, cs_advc_flags_s considers extended zones around buildings where
90!   first-order upwind scheme is applied for the horizontal advection terms,
91!   in order to overcome high concentration peaks due to stationary numerical
92!   oscillations caused by horizontal advection discretization.
93!
94! 4109 2019-07-22 17:00:34Z suehring
95! Slightly revise setting of boundary conditions at horizontal walls, use
96! data-structure offset index instead of pre-calculate it for each facing
97!
98! 4080 2019-07-09 18:17:37Z suehring
99! Restore accidantly removed limitation to positive values
100!
101! 4079 2019-07-09 18:04:41Z suehring
102! Application of monotonic flux limiter for the vertical scalar advection
103! up to the topography top (only for the cache-optimized version at the
104! moment).
105!
106! 4069 2019-07-01 14:05:51Z Giersch
107! Masked output running index mid has been introduced as a local variable to
108! avoid runtime error (Loop variable has been modified) in time_integration
109!
110! 4029 2019-06-14 14:04:35Z raasch
111! nest_chemistry option removed
112!
113! 4004 2019-05-24 11:32:38Z suehring
114! in subroutine chem_parin check emiss_lod / mod_emis only
115! when emissions_anthropogenic is activated in namelist (E.C. Chan)
116!
117! 3968 2019-05-13 11:04:01Z suehring
118! - added "emiss_lod" which serves the same function as "mode_emis"
119!   both will be synchronized with emiss_lod having pirority over
120!   mode_emis (see informational messages)
121! - modified existing error message and introduced new informational messages
122!    - CM0436 - now also applies to invalid LOD settings
123!    - CM0463 - emiss_lod take precedence in case of conflict with mod_emis
124!    - CM0464 - emiss_lod valued assigned based on mode_emis if undefined
125!
126! 3930 2019-04-24 14:57:18Z forkel
127! Changed chem_non_transport_physics to chem_non_advective_processes
128!
129!
130! 3929 2019-04-24 12:52:08Z banzhafs
131! Correct/complete module_interface introduction for chemistry model
132! Add subroutine chem_exchange_horiz_bounds
133! Bug fix deposition
134!
135! 3784 2019-03-05 14:16:20Z banzhafs
136! 2D output of emission fluxes
137!
138! 3784 2019-03-05 14:16:20Z banzhafs
139! Bugfix, uncomment erroneous commented variable used for dry deposition.
140! Bugfix in 3D emission output.
141!
142! 3784 2019-03-05 14:16:20Z banzhafs
143! Changes related to global restructuring of location messages and introduction
144! of additional debug messages
145!
146! 3784 2019-03-05 14:16:20Z banzhafs
147! some formatting of the deposition code
148!
149! 3784 2019-03-05 14:16:20Z banzhafs
150! some formatting
151!
152! 3784 2019-03-05 14:16:20Z banzhafs
153! added cs_mech to USE chem_gasphase_mod
154!
155! 3784 2019-03-05 14:16:20Z banzhafs
156! renamed get_mechanismname to get_mechanism_name
157! renamed do_emiss to emissions_anthropogenic and do_depo to deposition_dry (ecc)
158!
159! 3784 2019-03-05 14:16:20Z banzhafs
160! Unused variables removed/taken care of
161!
162!
163! 3784 2019-03-05 14:16:20Z forkel
164! Replaced READ from unit 10 by CALL get_mechanismname also in chem_header
165!
166!
167! 3783 2019-03-05 13:23:50Z forkel
168! Removed forgotte write statements an some unused variables (did not touch the
169! parts related to deposition)
170!
171!
172! 3780 2019-03-05 11:19:45Z forkel
173! Removed READ from unit 10, added CALL get_mechanismname
174!
175!
176! 3767 2019-02-27 08:18:02Z raasch
177! unused variable for file index removed from rrd-subroutines parameter list
178!
179! 3738 2019-02-12 17:00:45Z suehring
180! Clean-up debug prints
181!
182! 3737 2019-02-12 16:57:06Z suehring
183! Enable mesoscale offline nesting for chemistry variables as well as
184! initialization of chemistry via dynamic input file.
185!
186! 3719 2019-02-06 13:10:18Z kanani
187! Resolved cpu logpoint overlap with all progn.equations, moved cpu_log call
188! to prognostic_equations for better overview
189!
190! 3700 2019-01-26 17:03:42Z knoop
191! Some interface calls moved to module_interface + cleanup
192!
193! 3664 2019-01-09 14:00:37Z forkel
194! Replaced misplaced location message by @todo
195!
196!
197! 3654 2019-01-07 16:31:57Z suehring
198! Disable misplaced location message
199!
200! 3652 2019-01-07 15:29:59Z forkel
201! Checks added for chemistry mechanism, parameter chem_mechanism added (basit)
202!
203! 2718 2018-01-02 08:49:38Z maronga
204! Initial revision
205!
206!
207!
208!
209! Authors:
210! --------
211! @author Renate Forkel
212! @author Farah Kanani-Suehring
213! @author Klaus Ketelsen
214! @author Basit Khan
215! @author Sabine Banzhaf
216!
217!
218!------------------------------------------------------------------------------!
219! Description:
220! ------------
221!> Chemistry model for PALM-4U
222!> @todo Extend chem_species type by nspec and nvar as addititional elements (RF)
223!> @todo Check possibility to reduce dimension of chem_species%conc from nspec to nvar (RF)
224!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
225!>       allowed to use the chemistry model in a precursor run and additionally
226!>       not using it in a main run
227!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
228!> @todo Separate boundary conditions for each chem spcs to be implemented
229!> @todo Currently only total concentration are calculated. Resolved, parameterized
230!>       and chemistry fluxes although partially and some completely coded but
231!>       are not operational/activated in this version. bK.
232!> @todo slight differences in passive scalar and chem spcs when chem reactions
233!>       turned off. Need to be fixed. bK
234!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
235!> @todo chemistry error messages
236!
237!------------------------------------------------------------------------------!
238
239 MODULE chemistry_model_mod
240
241    USE advec_s_pw_mod,                                                                            &
242         ONLY:  advec_s_pw
243
244    USE advec_s_up_mod,                                                                            &
245         ONLY:  advec_s_up
246
247    USE advec_ws,                                                                                  &
248         ONLY:  advec_s_ws, ws_init_flags_scalar
249
250    USE diffusion_s_mod,                                                                           &
251         ONLY:  diffusion_s
252
253    USE kinds,                                                                                     &
254         ONLY:  iwp, wp
255
256    USE indices,                                                                                   &
257         ONLY:  advc_flags_s,                                                                      &
258                nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, nzb, nzt,            &
259                wall_flags_total_0
260
261    USE pegrid,                                                                                    &
262         ONLY: myid, threads_per_task
263
264    USE bulk_cloud_model_mod,                                                                      &
265         ONLY:  bulk_cloud_model
266
267    USE control_parameters,                                                                        &
268         ONLY:  bc_lr_cyc, bc_ns_cyc,                                                              &
269                bc_dirichlet_l,                                                                    &
270                bc_dirichlet_n,                                                                    &
271                bc_dirichlet_r,                                                                    &
272                bc_dirichlet_s,                                                                    &
273                bc_radiation_l,                                                                    &
274                bc_radiation_n,                                                                    &
275                bc_radiation_r,                                                                    &
276                bc_radiation_s,                                                                    &
277                debug_output,                                                                      &
278                dt_3d, humidity, initializing_actions, message_string,                             &
279                omega, tsc, intermediate_timestep_count, intermediate_timestep_count_max,          &
280                max_pr_user,                                                                       &
281                monotonic_limiter_z,                                                               &
282                scalar_advec,                                                                      &
283                timestep_scheme, use_prescribed_profile_data, ws_scheme_sca, air_chemistry
284
285    USE arrays_3d,                                                                                 &
286         ONLY:  exner, hyp, pt, q, ql, rdf_sc, tend, zu
287
288    USE chem_gasphase_mod,                                                                         &
289         ONLY:  atol, chem_gasphase_integrate, cs_mech, get_mechanism_name, nkppctrl,              &
290         nmaxfixsteps, nphot, nreact, nspec, nvar, phot_names, rtol, spc_names, t_steps, vl_dim
291
292    USE chem_modules
293
294    USE chem_photolysis_mod,                                                                       &
295        ONLY:  photolysis_control
296
297    USE cpulog,                                                                                    &
298        ONLY:  cpu_log, log_point_s
299
300    USE statistics
301
302    USE surface_mod,                                                                               &
303         ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
304
305    IMPLICIT NONE
306
307    PRIVATE
308    SAVE
309
310    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_1  !< pointer for swapping of timelevels for conc
311    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_2  !< pointer for swapping of timelevels for conc
312    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_3  !< pointer for swapping of timelevels for conc
313    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  spec_conc_av !< averaged concentrations of chemical species       
314    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET ::  freq_1       !< pointer for phtolysis frequncies
315                                                                            !< (only 1 timelevel required)
316    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl       !< 20 integer parameters for fine tuning KPP code
317                                                                            !< (e.g. solver type)
318    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl       !< 20 real parameters for fine tuning of KPP code
319                                                                            !< (e.g starting internal timestep of solver)
320!
321!-- Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
322!-- approproate for chemicals compounds since they may accumulate too much.
323!-- If no proper boundary conditions from a DYNAMIC input file are available,
324!-- de-cycling applies the initial profiles at the inflow boundaries for
325!-- Dirichlet boundary conditions
326    LOGICAL ::  decycle_chem_lr    = .FALSE.    !< switch for setting decycling in left-right direction
327    LOGICAL ::  decycle_chem_ns    = .FALSE.    !< switch for setting decycling in south-norht direction
328    CHARACTER (LEN=20), DIMENSION(4) ::  decycle_method = &
329         (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
330                              !< Decycling method at horizontal boundaries,
331                              !< 1=left, 2=right, 3=south, 4=north
332                              !< dirichlet = initial size distribution and
333                              !< chemical composition set for the ghost and
334                              !< first three layers
335                              !< neumann = zero gradient
336
337    REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
338
339!
340!-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
341    !
342    INTEGER(iwp), PARAMETER ::  nlu_dep = 15                   !< Number of DEPAC landuse classes (lu's)
343    INTEGER(iwp), PARAMETER ::  ncmp = 10                      !< Number of DEPAC gas components
344    INTEGER(iwp), PARAMETER ::  nposp = 69                     !< Number of possible species for deposition
345!
346!-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                             
347    INTEGER(iwp) ::  ilu_grass              = 1       
348    INTEGER(iwp) ::  ilu_arable             = 2       
349    INTEGER(iwp) ::  ilu_permanent_crops    = 3         
350    INTEGER(iwp) ::  ilu_coniferous_forest  = 4         
351    INTEGER(iwp) ::  ilu_deciduous_forest   = 5         
352    INTEGER(iwp) ::  ilu_water_sea          = 6       
353    INTEGER(iwp) ::  ilu_urban              = 7       
354    INTEGER(iwp) ::  ilu_other              = 8 
355    INTEGER(iwp) ::  ilu_desert             = 9 
356    INTEGER(iwp) ::  ilu_ice                = 10 
357    INTEGER(iwp) ::  ilu_savanna            = 11 
358    INTEGER(iwp) ::  ilu_tropical_forest    = 12 
359    INTEGER(iwp) ::  ilu_water_inland       = 13 
360    INTEGER(iwp) ::  ilu_mediterrean_scrub  = 14 
361    INTEGER(iwp) ::  ilu_semi_natural_veg   = 15 
362
363!
364!-- NH3/SO2 ratio regimes:
365    INTEGER(iwp), PARAMETER ::  iratns_low      = 1       !< low ratio NH3/SO2
366    INTEGER(iwp), PARAMETER ::  iratns_high     = 2       !< high ratio NH3/SO2
367    INTEGER(iwp), PARAMETER ::  iratns_very_low = 3       !< very low ratio NH3/SO2
368!
369!-- Default:
370    INTEGER, PARAMETER ::  iratns_default = iratns_low
371!
372!-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
373    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   =(/ 0.009, 0.009, 0.009, 0.006, 0.006, -999., -999., 0.009, -999.,      &
374         -999., 0.009, 0.006, -999., 0.009, 0.008/)*4.57
375!
376!-- Set temperatures per land use for f_temp
377    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999.,      &
378         12.0,  0.0, -999., 12.0,  8.0/)
379    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  topt = (/ 26.0, 26.0,  26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999.,      &
380         26.0, 20.0, -999., 26.0, 24.0 /)
381    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999.,      &
382         40.0, 35.0, -999., 40.0, 39.0 /)
383!
384!-- Set f_min:
385    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999., -999., 0.01, -999., -999., 0.01,  &
386         0.1, -999., 0.01, 0.04/)
387
388!
389!-- Set maximal conductance (m/s)
390!-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
391    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270., 300., 300., 140., 150., -999., -999., 270., -999., -999.,      &
392         270., 150., -999., 300., 422./)/41000.
393!
394!-- Set max, min for vapour pressure deficit vpd
395    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_max = (/1.3, 0.9, 0.9, 0.5, 1.0, -999., -999., 1.3, -999., -999., 1.3,      &
396         1.0, -999., 0.9, 2.8/) 
397    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_min = (/3.0, 2.8, 2.8, 3.0, 3.25, -999., -999., 3.0, -999., -999., 3.0,     &
398         3.25, -999., 2.8, 4.5/) 
399
400    PUBLIC nreact
401    PUBLIC nspec               !< number of gas phase chemical species including constant compound (e.g. N2)
402    PUBLIC nvar                !< number of variable gas phase chemical species (nvar <= nspec)
403    PUBLIC spc_names           !< names of gas phase chemical species (come from KPP) (come from KPP)
404    PUBLIC spec_conc_2 
405!   
406!-- Interface section
407    INTERFACE chem_3d_data_averaging
408       MODULE PROCEDURE chem_3d_data_averaging
409    END INTERFACE chem_3d_data_averaging
410
411    INTERFACE chem_boundary_conds
412       MODULE PROCEDURE chem_boundary_conds
413       MODULE PROCEDURE chem_boundary_conds_decycle
414    END INTERFACE chem_boundary_conds
415
416    INTERFACE chem_boundary_conditions
417       MODULE PROCEDURE chem_boundary_conditions
418    END INTERFACE chem_boundary_conditions
419
420    INTERFACE chem_check_data_output
421       MODULE PROCEDURE chem_check_data_output
422    END INTERFACE chem_check_data_output
423
424    INTERFACE chem_data_output_2d
425       MODULE PROCEDURE chem_data_output_2d
426    END INTERFACE chem_data_output_2d
427
428    INTERFACE chem_data_output_3d
429       MODULE PROCEDURE chem_data_output_3d
430    END INTERFACE chem_data_output_3d
431
432    INTERFACE chem_data_output_mask
433       MODULE PROCEDURE chem_data_output_mask
434    END INTERFACE chem_data_output_mask
435
436    INTERFACE chem_check_data_output_pr
437       MODULE PROCEDURE chem_check_data_output_pr
438    END INTERFACE chem_check_data_output_pr
439
440    INTERFACE chem_check_parameters
441       MODULE PROCEDURE chem_check_parameters
442    END INTERFACE chem_check_parameters
443
444    INTERFACE chem_define_netcdf_grid
445       MODULE PROCEDURE chem_define_netcdf_grid
446    END INTERFACE chem_define_netcdf_grid
447
448    INTERFACE chem_header
449       MODULE PROCEDURE chem_header
450    END INTERFACE chem_header
451
452    INTERFACE chem_init_arrays
453       MODULE PROCEDURE chem_init_arrays
454    END INTERFACE chem_init_arrays
455
456    INTERFACE chem_init
457       MODULE PROCEDURE chem_init
458    END INTERFACE chem_init
459
460    INTERFACE chem_init_profiles
461       MODULE PROCEDURE chem_init_profiles
462    END INTERFACE chem_init_profiles
463
464    INTERFACE chem_integrate
465       MODULE PROCEDURE chem_integrate_ij
466    END INTERFACE chem_integrate
467
468    INTERFACE chem_parin
469       MODULE PROCEDURE chem_parin
470    END INTERFACE chem_parin
471
472    INTERFACE chem_actions
473       MODULE PROCEDURE chem_actions
474       MODULE PROCEDURE chem_actions_ij
475    END INTERFACE chem_actions
476
477    INTERFACE chem_non_advective_processes
478       MODULE PROCEDURE chem_non_advective_processes
479       MODULE PROCEDURE chem_non_advective_processes_ij
480    END INTERFACE chem_non_advective_processes
481   
482    INTERFACE chem_exchange_horiz_bounds
483       MODULE PROCEDURE chem_exchange_horiz_bounds
484    END INTERFACE chem_exchange_horiz_bounds   
485
486    INTERFACE chem_prognostic_equations
487       MODULE PROCEDURE chem_prognostic_equations
488       MODULE PROCEDURE chem_prognostic_equations_ij
489    END INTERFACE chem_prognostic_equations
490
491    INTERFACE chem_rrd_local
492       MODULE PROCEDURE chem_rrd_local
493    END INTERFACE chem_rrd_local
494
495    INTERFACE chem_statistics
496       MODULE PROCEDURE chem_statistics
497    END INTERFACE chem_statistics
498
499    INTERFACE chem_swap_timelevel
500       MODULE PROCEDURE chem_swap_timelevel
501    END INTERFACE chem_swap_timelevel
502
503    INTERFACE chem_wrd_local
504       MODULE PROCEDURE chem_wrd_local
505    END INTERFACE chem_wrd_local
506
507    INTERFACE chem_depo
508       MODULE PROCEDURE chem_depo
509    END INTERFACE chem_depo
510
511    INTERFACE drydepos_gas_depac
512       MODULE PROCEDURE drydepos_gas_depac
513    END INTERFACE drydepos_gas_depac
514
515    INTERFACE rc_special
516       MODULE PROCEDURE rc_special
517    END INTERFACE rc_special
518
519    INTERFACE  rc_gw
520       MODULE PROCEDURE rc_gw
521    END INTERFACE rc_gw
522
523    INTERFACE rw_so2
524       MODULE PROCEDURE rw_so2 
525    END INTERFACE rw_so2
526
527    INTERFACE rw_nh3_sutton
528       MODULE PROCEDURE rw_nh3_sutton
529    END INTERFACE rw_nh3_sutton
530
531    INTERFACE rw_constant
532       MODULE PROCEDURE rw_constant
533    END INTERFACE rw_constant
534
535    INTERFACE rc_gstom
536       MODULE PROCEDURE rc_gstom
537    END INTERFACE rc_gstom
538
539    INTERFACE rc_gstom_emb
540       MODULE PROCEDURE rc_gstom_emb
541    END INTERFACE rc_gstom_emb
542
543    INTERFACE par_dir_diff
544       MODULE PROCEDURE par_dir_diff
545    END INTERFACE par_dir_diff
546
547    INTERFACE rc_get_vpd
548       MODULE PROCEDURE rc_get_vpd
549    END INTERFACE rc_get_vpd
550
551    INTERFACE rc_gsoil_eff
552       MODULE PROCEDURE rc_gsoil_eff
553    END INTERFACE rc_gsoil_eff
554
555    INTERFACE rc_rinc
556       MODULE PROCEDURE rc_rinc
557    END INTERFACE rc_rinc
558
559    INTERFACE rc_rctot
560       MODULE PROCEDURE rc_rctot 
561    END INTERFACE rc_rctot
562
563!    INTERFACE rc_comp_point_rc_eff
564!       MODULE PROCEDURE rc_comp_point_rc_eff
565!    END INTERFACE rc_comp_point_rc_eff
566
567    INTERFACE drydepo_aero_zhang_vd
568       MODULE PROCEDURE drydepo_aero_zhang_vd
569    END INTERFACE drydepo_aero_zhang_vd
570
571    INTERFACE get_rb_cell
572       MODULE PROCEDURE  get_rb_cell
573    END INTERFACE get_rb_cell
574
575
576
577    PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_boundary_conditions, &
578            chem_boundary_conds_decycle, chem_check_data_output,              &
579         chem_check_data_output_pr, chem_check_parameters,                    &
580         chem_data_output_2d, chem_data_output_3d, chem_data_output_mask,     &
581         chem_define_netcdf_grid, chem_header, chem_init, chem_init_arrays,   &
582         chem_init_profiles, chem_integrate, chem_parin,                      &
583         chem_actions, chem_prognostic_equations, chem_rrd_local,             &
584         chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo,     &
585         chem_non_advective_processes, chem_exchange_horiz_bounds
586
587 CONTAINS
588
589
590!------------------------------------------------------------------------------!
591! Description:
592! ------------
593!> Subroutine for averaging 3D data of chemical species. Due to the fact that
594!> the averaged chem arrays are allocated in chem_init, no if-query concerning
595!> the allocation is required (in any mode). Attention: If you just specify an
596!> averaged output quantity in the _p3dr file during restarts the first output
597!> includes the time between the beginning of the restart run and the first
598!> output time (not necessarily the whole averaging_interval you have
599!> specified in your _p3d/_p3dr file )
600!------------------------------------------------------------------------------!
601 SUBROUTINE chem_3d_data_averaging( mode, variable )
602
603
604    USE control_parameters
605
606    CHARACTER (LEN=*) ::  mode     !<
607    CHARACTER (LEN=*) ::  variable !<
608
609    LOGICAL ::  match_def  !< flag indicating default-type surface
610    LOGICAL ::  match_lsm  !< flag indicating natural-type surface
611    LOGICAL ::  match_usm  !< flag indicating urban-type surface
612
613    INTEGER(iwp) ::  i                  !< grid index x direction
614    INTEGER(iwp) ::  j                  !< grid index y direction
615    INTEGER(iwp) ::  k                  !< grid index z direction
616    INTEGER(iwp) ::  m                  !< running index surface type
617    INTEGER(iwp) ::  lsp               !< running index for chem spcs
618
619    IF ( (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_')  )  THEN
620
621       IF ( mode == 'allocate' )  THEN
622
623          DO  lsp = 1, nspec
624             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
625                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
626                chem_species(lsp)%conc_av = 0.0_wp
627             ENDIF
628          ENDDO
629
630       ELSEIF ( mode == 'sum' )  THEN
631
632          DO  lsp = 1, nspec
633             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
634                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
635                DO  i = nxlg, nxrg
636                   DO  j = nysg, nyng
637                      DO  k = nzb, nzt+1
638                         chem_species(lsp)%conc_av(k,j,i) =                              &
639                                           chem_species(lsp)%conc_av(k,j,i) +            &
640                                           chem_species(lsp)%conc(k,j,i)
641                      ENDDO
642                   ENDDO
643                ENDDO
644             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
645                DO  i = nxl, nxr
646                   DO  j = nys, nyn
647                      match_def = surf_def_h(0)%start_index(j,i) <=                      &
648                           surf_def_h(0)%end_index(j,i)
649                      match_lsm = surf_lsm_h%start_index(j,i) <=                         &
650                           surf_lsm_h%end_index(j,i)
651                      match_usm = surf_usm_h%start_index(j,i) <=                         &
652                           surf_usm_h%end_index(j,i)
653
654                      IF ( match_def )  THEN
655                         m = surf_def_h(0)%end_index(j,i)
656                         chem_species(lsp)%cssws_av(j,i) =                               &
657                              chem_species(lsp)%cssws_av(j,i) + &
658                              surf_def_h(0)%cssws(lsp,m)
659                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
660                         m = surf_lsm_h%end_index(j,i)
661                         chem_species(lsp)%cssws_av(j,i) =                               &
662                              chem_species(lsp)%cssws_av(j,i) + &
663                              surf_lsm_h%cssws(lsp,m)
664                      ELSEIF ( match_usm )  THEN
665                         m = surf_usm_h%end_index(j,i)
666                         chem_species(lsp)%cssws_av(j,i) =                               &
667                              chem_species(lsp)%cssws_av(j,i) + &
668                              surf_usm_h%cssws(lsp,m)
669                      ENDIF
670                   ENDDO
671                ENDDO
672             ENDIF
673          ENDDO
674
675       ELSEIF ( mode == 'average' )  THEN
676
677          DO  lsp = 1, nspec
678             IF ( TRIM( variable(1:3) ) == 'kc_' .AND. &
679                  TRIM( variable(4:) ) == TRIM( chem_species(lsp)%name ) )  THEN
680                DO  i = nxlg, nxrg
681                   DO  j = nysg, nyng
682                      DO  k = nzb, nzt+1
683                         chem_species(lsp)%conc_av(k,j,i) =                              &
684                             chem_species(lsp)%conc_av(k,j,i) /                          &
685                             REAL( average_count_3d, KIND=wp )
686                      ENDDO
687                   ENDDO
688                ENDDO
689
690             ELSEIF ( TRIM( variable(4:) ) == TRIM( 'cssws*' ) )  THEN
691                DO  i = nxlg, nxrg
692                   DO  j = nysg, nyng
693                      chem_species(lsp)%cssws_av(j,i) =                                  &
694                      chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
695                   ENDDO
696                ENDDO
697                CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
698             ENDIF
699          ENDDO
700       ENDIF
701
702    ENDIF
703
704 END SUBROUTINE chem_3d_data_averaging
705
706   
707!------------------------------------------------------------------------------!
708! Description:
709! ------------
710!> Subroutine to initialize and set all boundary conditions for chemical species
711!------------------------------------------------------------------------------!
712 SUBROUTINE chem_boundary_conds( mode )                                           
713
714    USE control_parameters,                                                    & 
715        ONLY:  bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s
716
717    USE arrays_3d,                                                             &     
718        ONLY:  dzu                                               
719
720    USE surface_mod,                                                           &
721        ONLY:  bc_h                                                           
722
723    CHARACTER (LEN=*), INTENT(IN) ::  mode
724    INTEGER(iwp) ::  i                            !< grid index x direction.
725    INTEGER(iwp) ::  j                            !< grid index y direction.
726    INTEGER(iwp) ::  k                            !< grid index z direction.
727    INTEGER(iwp) ::  l                            !< running index boundary type, for up- and downward-facing walls.
728    INTEGER(iwp) ::  m                            !< running index surface elements.
729    INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
730
731
732    SELECT CASE ( TRIM( mode ) )       
733       CASE ( 'init' )
734
735          IF ( bc_cs_b == 'dirichlet' )  THEN
736             ibc_cs_b = 0 
737          ELSEIF ( bc_cs_b == 'neumann' )  THEN
738             ibc_cs_b = 1 
739          ELSE
740             message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 
741             CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )
742          ENDIF                                                                 
743!
744!--       Set Integer flags and check for possible erroneous settings for top
745!--       boundary condition.
746          IF ( bc_cs_t == 'dirichlet' )  THEN
747             ibc_cs_t = 0 
748          ELSEIF ( bc_cs_t == 'neumann' )  THEN
749             ibc_cs_t = 1
750          ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
751             ibc_cs_t = 2
752          ELSEIF ( bc_cs_t == 'nested' )  THEN         
753             ibc_cs_t = 3
754          ELSE
755             message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'     
756             CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
757          ENDIF
758
759       CASE ( 'set_bc_bottomtop' )                   
760!
761!--       Boundary condtions for chemical species at horizontal walls     
762          DO  lsp = 1, nspec                                                     
763             IF ( ibc_cs_b == 0 )  THEN
764                DO  l = 0, 1 
765                   !$OMP PARALLEL DO PRIVATE( i, j, k )
766                   DO  m = 1, bc_h(l)%ns
767                       i = bc_h(l)%i(m)           
768                       j = bc_h(l)%j(m)
769                       k = bc_h(l)%k(m)
770                      chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
771                                      chem_species(lsp)%conc(k+bc_h(l)%koff,j,i) 
772                   ENDDO                                       
773                ENDDO                                       
774
775             ELSEIF ( ibc_cs_b == 1 )  THEN
776!
777!--             in boundary_conds there is som extra loop over m here for passive tracer
778                DO  l = 0, 1
779                   !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
780                   DO m = 1, bc_h(l)%ns
781                      i = bc_h(l)%i(m)           
782                      j = bc_h(l)%j(m)
783                      k = bc_h(l)%k(m)
784                      chem_species(lsp)%conc_p(k+bc_h(l)%koff,j,i) =           &
785                                         chem_species(lsp)%conc_p(k,j,i)
786
787                   ENDDO
788                ENDDO
789             ENDIF
790       ENDDO    ! end lsp loop 
791!
792!--    Top boundary conditions for chemical species - Should this not be done for all species?
793          IF ( ibc_cs_t == 0 )  THEN                   
794             DO  lsp = 1, nspec
795                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)       
796             ENDDO
797          ELSEIF ( ibc_cs_t == 1 )  THEN
798             DO  lsp = 1, nspec
799                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
800             ENDDO
801          ELSEIF ( ibc_cs_t == 2 )  THEN
802             DO  lsp = 1, nspec
803                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
804             ENDDO
805          ENDIF
806
807       CASE ( 'set_bc_lateral' )                       
808!
809!--             Lateral boundary conditions for chem species at inflow boundary
810!--             are automatically set when chem_species concentration is
811!--             initialized. The initially set value at the inflow boundary is not
812!--             touched during time integration, hence, this boundary value remains
813!--             at a constant value, which is the concentration that flows into the
814!--             domain.                                                           
815!--             Lateral boundary conditions for chem species at outflow boundary
816
817          IF ( bc_radiation_s )  THEN
818             DO  lsp = 1, nspec
819                chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
820             ENDDO
821          ELSEIF ( bc_radiation_n )  THEN
822             DO  lsp = 1, nspec
823                chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
824             ENDDO
825          ELSEIF ( bc_radiation_l )  THEN
826             DO  lsp = 1, nspec
827                chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
828             ENDDO
829          ELSEIF ( bc_radiation_r )  THEN
830             DO  lsp = 1, nspec
831                chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
832             ENDDO
833          ENDIF
834
835    END SELECT
836
837 END SUBROUTINE chem_boundary_conds
838
839!------------------------------------------------------------------------------!
840! Description:
841! ------------
842!> Subroutine for boundary conditions
843!------------------------------------------------------------------------------!
844 SUBROUTINE chem_boundary_conditions
845
846    IMPLICIT NONE
847
848    INTEGER(iwp) ::  lsp             !<
849    INTEGER(iwp) ::  lsp_usr         !<
850
851!
852!-- Top/bottom boundary conditions for chemical species
853    CALL chem_boundary_conds( 'set_bc_bottomtop' )
854!
855!-- Lateral boundary conditions for chemical species
856    CALL chem_boundary_conds( 'set_bc_lateral' )
857
858!
859!--  Boundary conditions for prognostic quantitites of other modules:
860!--  Here, only decycling is carried out
861     DO  lsp = 1, nvar
862        lsp_usr = 1
863        DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
864           IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
865              CALL chem_boundary_conds( chem_species(lsp)%conc_p,                          &
866                                        chem_species(lsp)%conc_pr_init )
867           ENDIF
868           lsp_usr = lsp_usr + 1
869        ENDDO
870     ENDDO
871
872
873 END SUBROUTINE chem_boundary_conditions
874
875!------------------------------------------------------------------------------!
876! Description:
877! ------------
878!> Boundary conditions for prognostic variables in chemistry: decycling in the
879!> x-direction-
880!> Decycling of chemistry variables: Dirichlet BCs with cyclic is frequently not
881!> approproate for chemicals compounds since they may accumulate too much.
882!> If no proper boundary conditions from a DYNAMIC input file are available,
883!> de-cycling applies the initial profiles at the inflow boundaries for
884!> Dirichlet boundary conditions
885!------------------------------------------------------------------------------!
886 SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init )
887
888    USE control_parameters,                                                                         &
889        ONLY:  nesting_offline
890
891    INTEGER(iwp) ::  boundary  !<
892    INTEGER(iwp) ::  ee        !<
893    INTEGER(iwp) ::  copied    !<
894    INTEGER(iwp) ::  i         !<
895    INTEGER(iwp) ::  j         !<
896    INTEGER(iwp) ::  k         !<
897    INTEGER(iwp) ::  ss        !<
898
899    REAL(wp), DIMENSION(nzb:nzt+1) ::  cs_pr_init
900    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  cs_3d
901    REAL(wp) ::  flag          !< flag to mask topography grid points
902
903
904    flag = 0.0_wp
905!
906!-- Skip input if forcing from a larger-scale model is applied
907    IF ( nesting_offline  .AND.  nesting_offline_chem )  RETURN
908!
909!-- Left and right boundaries
910    IF ( decycle_chem_lr  .AND.  bc_lr_cyc )  THEN
911
912       DO  boundary = 1, 2
913
914          IF ( decycle_method(boundary) == 'dirichlet' )  THEN
915!
916!--          Initial profile is copied to ghost and first three layers
917             ss = 1
918             ee = 0
919             IF ( boundary == 1  .AND.  nxl == 0 )  THEN
920                ss = nxlg
921                ee = nxl-1
922             ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
923                ss = nxr+1
924                ee = nxrg
925             ENDIF
926
927             DO  i = ss, ee
928                DO  j = nysg, nyng
929                   DO  k = nzb+1, nzt
930                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
931                      cs_3d(k,j,i) = cs_pr_init(k) * flag
932                   ENDDO
933                ENDDO
934             ENDDO
935
936          ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
937!
938!--          The value at the boundary is copied to the ghost layers to simulate
939!--          an outlet with zero gradient
940             ss = 1
941             ee = 0
942             IF ( boundary == 1  .AND.  nxl == 0 )  THEN
943                ss = nxlg
944                ee = nxl-1
945                copied = nxl
946             ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
947                ss = nxr+1
948                ee = nxrg
949                copied = nxr
950             ENDIF
951
952             DO  i = ss, ee
953                DO  j = nysg, nyng
954                   DO  k = nzb+1, nzt
955                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
956                      cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
957                   ENDDO
958                ENDDO
959             ENDDO
960
961          ELSE
962             WRITE(message_string,*)                                           &
963                  'unknown decycling method: decycle_method (', &
964                  boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
965             CALL message( 'chem_boundary_conds_decycle', 'CM0431',           &
966                  1, 2, 0, 6, 0 )
967          ENDIF
968       ENDDO
969    ENDIF
970!
971!-- South and north boundaries
972    IF ( decycle_chem_ns  .AND.  bc_ns_cyc )  THEN
973
974       DO  boundary = 3, 4
975
976          IF ( decycle_method(boundary) == 'dirichlet' )  THEN
977!
978!--          Initial profile is copied to ghost and first three layers
979             ss = 1
980             ee = 0
981             IF ( boundary == 3  .AND.  nys == 0 )  THEN
982                ss = nysg
983                ee = nys-1
984             ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
985                ss = nyn+1
986                ee = nyng
987             ENDIF
988
989             DO  i = nxlg, nxrg
990                DO  j = ss, ee
991                   DO  k = nzb+1, nzt
992                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
993                      cs_3d(k,j,i) = cs_pr_init(k) * flag
994                   ENDDO
995                ENDDO
996             ENDDO
997
998
999          ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
1000!
1001!--          The value at the boundary is copied to the ghost layers to simulate
1002!--          an outlet with zero gradient
1003             ss = 1
1004             ee = 0
1005             IF ( boundary == 3  .AND.  nys == 0 )  THEN
1006                ss = nysg
1007                ee = nys-1
1008                copied = nys
1009             ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
1010                ss = nyn+1
1011                ee = nyng
1012                copied = nyn
1013             ENDIF
1014
1015             DO  i = nxlg, nxrg
1016                DO  j = ss, ee
1017                   DO  k = nzb+1, nzt
1018                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
1019                      cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
1020                   ENDDO
1021                ENDDO
1022             ENDDO
1023
1024          ELSE
1025             WRITE(message_string,*)                                           &
1026                  'unknown decycling method: decycle_method (', &
1027                  boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
1028             CALL message( 'chem_boundary_conds_decycle', 'CM0432',           &
1029                  1, 2, 0, 6, 0 )
1030          ENDIF
1031       ENDDO
1032    ENDIF
1033
1034
1035 END SUBROUTINE chem_boundary_conds_decycle
1036
1037
1038!------------------------------------------------------------------------------!
1039! Description:
1040! ------------
1041!> Subroutine for checking data output for chemical species
1042!------------------------------------------------------------------------------!
1043 SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
1044
1045
1046    CHARACTER (LEN=*) ::  unit     !<
1047    CHARACTER (LEN=*) ::  var      !<
1048
1049    INTEGER(iwp) ::  i
1050    INTEGER(iwp) ::  lsp
1051    INTEGER(iwp) ::  ilen
1052    INTEGER(iwp) ::  k
1053
1054    CHARACTER(LEN=16)    ::  spec_name
1055
1056!
1057!-- Next statement is to avoid compiler warnings about unused variables
1058    IF ( ( i + ilen + k ) > 0  .OR.  var(1:1) == ' ' )  CONTINUE
1059
1060    unit = 'illegal'
1061
1062    spec_name = TRIM( var(4:) )             !< var 1:3 is 'kc_' or 'em_'.
1063
1064    IF ( TRIM( var(1:3) ) == 'em_' )  THEN
1065       DO  lsp=1,nspec
1066          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1067             unit = 'mol m-2 s-1'
1068          ENDIF
1069!
1070!--       It is possible to plant PM10 and PM25 into the gasphase chemistry code
1071!--       as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
1072!--       set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
1073          IF (spec_name(1:2) == 'PM')  THEN
1074             unit = 'kg m-2 s-1'
1075          ENDIF
1076       ENDDO
1077
1078    ELSE
1079
1080       DO  lsp=1,nspec
1081          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1082             unit = 'ppm'
1083          ENDIF
1084!
1085!--            It is possible  to plant PM10 and PM25 into the gasphase chemistry code
1086!--            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
1087!--            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
1088          IF (spec_name(1:2) == 'PM')  THEN 
1089            unit = 'kg m-3'
1090          ENDIF
1091       ENDDO
1092
1093       DO  lsp=1,nphot
1094          IF (TRIM( spec_name ) == TRIM( phot_frequen(lsp)%name ) )  THEN
1095             unit = 'sec-1'
1096          ENDIF
1097       ENDDO
1098    ENDIF
1099
1100
1101    RETURN
1102 END SUBROUTINE chem_check_data_output
1103
1104
1105!------------------------------------------------------------------------------!
1106! Description:
1107! ------------
1108!> Subroutine for checking data output of profiles for chemistry model
1109!------------------------------------------------------------------------------!
1110 SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
1111
1112    USE arrays_3d
1113
1114    USE control_parameters,                                                    &
1115        ONLY:  data_output_pr, message_string
1116
1117    USE profil_parameter
1118
1119    USE statistics
1120
1121
1122    CHARACTER (LEN=*) ::  unit     !<
1123    CHARACTER (LEN=*) ::  variable !<
1124    CHARACTER (LEN=*) ::  dopr_unit
1125    CHARACTER (LEN=16) ::  spec_name
1126
1127    INTEGER(iwp) ::  var_count, lsp  !<
1128
1129
1130    spec_name = TRIM( variable(4:) )   
1131
1132       IF (  .NOT.  air_chemistry )  THEN
1133          message_string = 'data_output_pr = ' //                        &
1134          TRIM( data_output_pr(var_count) ) // ' is not imp' // &
1135                      'lemented for air_chemistry = .FALSE.'
1136          CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )             
1137
1138       ELSE
1139          DO  lsp = 1, nspec
1140             IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) )  THEN
1141                cs_pr_count = cs_pr_count+1
1142                cs_pr_index(cs_pr_count) = lsp
1143                dopr_index(var_count) = pr_palm+cs_pr_count 
1144                dopr_unit  = 'ppm'
1145                IF (spec_name(1:2) == 'PM')  THEN
1146                     dopr_unit = 'kg m-3'
1147                ENDIF
1148                   hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
1149                   unit = dopr_unit 
1150             ENDIF
1151          ENDDO
1152       ENDIF
1153
1154 END SUBROUTINE chem_check_data_output_pr
1155
1156
1157!------------------------------------------------------------------------------!
1158! Description:
1159! ------------
1160!> Check parameters routine for chemistry_model_mod
1161!------------------------------------------------------------------------------!
1162 SUBROUTINE chem_check_parameters
1163
1164
1165    LOGICAL  ::  found
1166    INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs
1167    INTEGER (iwp) ::  lsp          !< running index for chem spcs.
1168!
1169!-- check for chemical reactions status
1170    IF ( chem_gasphase_on )  THEN
1171       message_string = 'Chemical reactions: ON'
1172       CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
1173    ELSEIF ( .NOT. (chem_gasphase_on) )  THEN
1174       message_string = 'Chemical reactions: OFF'
1175       CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
1176    ENDIF
1177!
1178!-- check for chemistry time-step
1179    IF ( call_chem_at_all_substeps )  THEN
1180       message_string = 'Chemistry is calculated at all meteorology time-step'
1181       CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
1182    ELSEIF ( .not. (call_chem_at_all_substeps) )  THEN
1183       message_string = 'Sub-time-steps are skipped for chemistry time-steps'
1184       CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
1185    ENDIF
1186!
1187!-- check for photolysis scheme
1188    IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
1189       message_string = 'Incorrect photolysis scheme selected, please check spelling'
1190       CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
1191    ENDIF
1192!
1193!-- check for  decycling of chem species
1194    IF ( (.NOT. any(decycle_method == 'neumann') ) .AND. (.NOT. any(decycle_method == 'dirichlet') ) )  THEN
1195       message_string = 'Incorrect boundary conditions. Only neumann or '   &
1196                // 'dirichlet &available for decycling chemical species '
1197       CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
1198    ENDIF
1199!
1200!-- check for chemical mechanism used
1201    CALL get_mechanism_name
1202    IF ( chem_mechanism /= TRIM( cs_mech ) )  THEN
1203       message_string = 'Incorrect chemistry mechanism selected, check spelling in namelist and/or chem_gasphase_mod'
1204       CALL message( 'chem_check_parameters', 'CM0462', 1, 2, 0, 6, 0 )
1205    ENDIF
1206!
1207!-- If nesting_chem = .F., set top boundary condition to its default value
1208    IF ( .NOT. nesting_chem  .AND.  ibc_cs_t == 3  )  THEN
1209       ibc_cs_t = 2
1210       bc_cs_t = 'initial_gradient'
1211    ENDIF
1212!
1213!-- chem_check_parameters is called before the array chem_species is allocated!
1214!-- temporary switch of this part of the check
1215!    RETURN                !bK commented
1216
1217    CALL chem_init_internal
1218!
1219!-- check for initial chem species input
1220    lsp_usr = 1
1221    lsp     = 1
1222    DO WHILE ( cs_name (lsp_usr) /= 'novalue')
1223       found = .FALSE.
1224       DO  lsp = 1, nvar
1225          IF ( TRIM( cs_name (lsp_usr) ) == TRIM( chem_species(lsp)%name) )  THEN
1226             found = .TRUE.
1227             EXIT
1228          ENDIF
1229       ENDDO
1230       IF ( .NOT.  found )  THEN
1231          message_string = 'Unused/incorrect input for initial surface value: ' //     &
1232                            TRIM( cs_name(lsp_usr) )
1233          CALL message( 'chem_check_parameters', 'CM0427', 1, 2, 0, 6, 0 )
1234       ENDIF
1235       lsp_usr = lsp_usr + 1
1236    ENDDO
1237!
1238!-- check for surface  emission flux chem species
1239    lsp_usr = 1
1240    lsp     = 1
1241    DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
1242       found = .FALSE.
1243       DO  lsp = 1, nvar
1244          IF ( TRIM( surface_csflux_name (lsp_usr) ) == TRIM( chem_species(lsp)%name ) )  THEN
1245             found = .TRUE.
1246             EXIT
1247          ENDIF
1248       ENDDO
1249       IF ( .NOT.  found )  THEN
1250          message_string = 'Unused/incorrect input of chemical species for surface emission fluxes: '  &
1251                            // TRIM( surface_csflux_name(lsp_usr) )
1252          CALL message( 'chem_check_parameters', 'CM0428', 1, 2, 0, 6, 0 )
1253       ENDIF
1254       lsp_usr = lsp_usr + 1
1255    ENDDO
1256
1257 END SUBROUTINE chem_check_parameters
1258
1259
1260!------------------------------------------------------------------------------!
1261! Description:
1262! ------------
1263!> Subroutine defining 2D output variables for chemical species
1264!> @todo: Remove "mode" from argument list, not used.
1265!------------------------------------------------------------------------------!
1266 SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf,   &
1267                                  two_d, nzb_do, nzt_do, fill_value )
1268
1269
1270    CHARACTER (LEN=*) ::  grid       !<
1271    CHARACTER (LEN=*) ::  mode       !<
1272    CHARACTER (LEN=*) ::  variable   !<
1273    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
1274    INTEGER(iwp) ::  nzb_do          !< lower limit of the domain (usually nzb)
1275    INTEGER(iwp) ::  nzt_do          !< upper limit of the domain (usually nzt+1)
1276    LOGICAL      ::  found           !<
1277    LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross sections)
1278    REAL(wp)     ::  fill_value
1279    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf 
1280
1281!
1282!-- local variables.
1283    CHARACTER(LEN=16)    ::  spec_name
1284    INTEGER(iwp) ::  lsp
1285    INTEGER(iwp) ::  i               !< grid index along x-direction
1286    INTEGER(iwp) ::  j               !< grid index along y-direction
1287    INTEGER(iwp) ::  k               !< grid index along z-direction
1288    INTEGER(iwp) ::  m               !< running indices for surfaces
1289    INTEGER(iwp) ::  char_len        !< length of a character string
1290!
1291!-- Next statement is to avoid compiler warnings about unused variables
1292    IF ( mode(1:1) == ' '  .OR.  two_d )  CONTINUE
1293
1294    found = .FALSE.
1295    char_len  = LEN_TRIM( variable )
1296
1297    spec_name = TRIM( variable(4:char_len-3) )
1298!
1299!-- Output of emission values, i.e. surface fluxes cssws.
1300    IF ( variable(1:3) == 'em_' )  THEN
1301
1302       local_pf = 0.0_wp
1303
1304       DO  lsp = 1, nvar
1305          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1306!
1307!--          No average output for now.
1308             DO  m = 1, surf_lsm_h%ns
1309                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1) =              &
1310                              local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),nzb+1)  &
1311                                            + surf_lsm_h%cssws(lsp,m)
1312             ENDDO
1313             DO  m = 1, surf_usm_h%ns
1314                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1) =           &
1315                              local_pf(surf_usm_h%i(m),surf_usm_h%j(m),nzb+1)  &
1316                                            + surf_usm_h%cssws(lsp,m)
1317             ENDDO
1318             grid = 'zu'
1319             found = .TRUE.
1320          ENDIF
1321       ENDDO
1322
1323    ELSE
1324
1325       DO  lsp=1,nspec
1326          IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name )  .AND.                           &
1327                ( (variable(char_len-2:) == '_xy')  .OR.                                           &
1328                  (variable(char_len-2:) == '_xz')  .OR.                                           &
1329                  (variable(char_len-2:) == '_yz') ) )  THEN             
1330!
1331!--   todo: remove or replace by "CALL message" mechanism (kanani)
1332!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM( variable )  //       &
1333!                                                             TRIM( chem_species(lsp)%name )       
1334             IF (av == 0)  THEN
1335                DO  i = nxl, nxr
1336                   DO  j = nys, nyn
1337                      DO  k = nzb_do, nzt_do
1338                           local_pf(i,j,k) = MERGE(                                                &
1339                                              chem_species(lsp)%conc(k,j,i),                       &
1340                                              REAL( fill_value, KIND = wp ),                       &
1341                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
1342                      ENDDO
1343                   ENDDO
1344                ENDDO
1345
1346             ELSE
1347                DO  i = nxl, nxr
1348                   DO  j = nys, nyn
1349                      DO  k = nzb_do, nzt_do
1350                           local_pf(i,j,k) = MERGE(                                                &
1351                                              chem_species(lsp)%conc_av(k,j,i),                    &
1352                                              REAL( fill_value, KIND = wp ),                       &
1353                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
1354                      ENDDO
1355                   ENDDO
1356                ENDDO
1357             ENDIF
1358             grid = 'zu'           
1359             found = .TRUE.
1360          ENDIF
1361       ENDDO
1362    ENDIF
1363
1364    RETURN
1365
1366 END SUBROUTINE chem_data_output_2d     
1367
1368
1369!------------------------------------------------------------------------------!
1370! Description:
1371! ------------
1372!> Subroutine defining 3D output variables for chemical species
1373!------------------------------------------------------------------------------!
1374 SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1375
1376
1377    USE surface_mod
1378
1379    CHARACTER (LEN=*)    ::  variable     !<
1380    INTEGER(iwp)         ::  av           !<
1381    INTEGER(iwp) ::  nzb_do               !< lower limit of the data output (usually 0)
1382    INTEGER(iwp) ::  nzt_do               !< vertical upper limit of the data output (usually nz_do3d)
1383
1384    LOGICAL      ::  found                !<
1385
1386    REAL(wp)             ::  fill_value   !<
1387    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf
1388!
1389!-- local variables
1390    CHARACTER(LEN=16)    ::  spec_name
1391    INTEGER(iwp)         ::  i
1392    INTEGER(iwp)         ::  j
1393    INTEGER(iwp)         ::  k
1394    INTEGER(iwp)         ::  m       !< running indices for surfaces
1395    INTEGER(iwp)         ::  l
1396    INTEGER(iwp)         ::  lsp     !< running index for chem spcs
1397
1398
1399    found = .FALSE.
1400    IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) )  THEN
1401       RETURN
1402    ENDIF
1403
1404    spec_name = TRIM( variable(4:) )
1405
1406    IF ( variable(1:3) == 'em_' )  THEN
1407
1408       DO  lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
1409          IF ( TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1410         
1411             local_pf = 0.0_wp
1412!
1413!--          no average for now
1414             DO  m = 1, surf_usm_h%ns
1415                local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
1416                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
1417             ENDDO
1418             DO  m = 1, surf_lsm_h%ns
1419                local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
1420                  local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
1421             ENDDO
1422             DO  l = 0, 3
1423                DO  m = 1, surf_usm_v(l)%ns
1424                   local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
1425                     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)
1426                ENDDO
1427                DO  m = 1, surf_lsm_v(l)%ns
1428                   local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
1429                      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)
1430                ENDDO
1431             ENDDO
1432             found = .TRUE.
1433          ENDIF
1434       ENDDO
1435    ELSE
1436      DO  lsp = 1, nspec
1437         IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN
1438!
1439!--   todo: remove or replace by "CALL message" mechanism (kanani)
1440!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
1441!                                                           TRIM( chem_species(lsp)%name )       
1442            IF (av == 0)  THEN
1443               DO  i = nxl, nxr
1444                  DO  j = nys, nyn
1445                     DO  k = nzb_do, nzt_do
1446                         local_pf(i,j,k) = MERGE(                             &
1447                                             chem_species(lsp)%conc(k,j,i),   &
1448                                             REAL( fill_value, KIND = wp ),   &
1449                                             BTEST( wall_flags_total_0(k,j,i), 0 ) )
1450                     ENDDO
1451                  ENDDO
1452               ENDDO
1453
1454            ELSE
1455
1456               DO  i = nxl, nxr
1457                  DO  j = nys, nyn
1458                     DO  k = nzb_do, nzt_do
1459                         local_pf(i,j,k) = MERGE(                             &
1460                                             chem_species(lsp)%conc_av(k,j,i),&
1461                                             REAL( fill_value, KIND = wp ),   &
1462                                             BTEST( wall_flags_total_0(k,j,i), 0 ) )
1463                     ENDDO
1464                  ENDDO
1465               ENDDO
1466            ENDIF
1467            found = .TRUE.
1468         ENDIF
1469      ENDDO
1470    ENDIF
1471
1472    RETURN
1473
1474 END SUBROUTINE chem_data_output_3d
1475
1476
1477!------------------------------------------------------------------------------!
1478! Description:
1479! ------------
1480!> Subroutine defining mask output variables for chemical species
1481!------------------------------------------------------------------------------!
1482 SUBROUTINE chem_data_output_mask( av, variable, found, local_pf, mid )
1483
1484
1485    USE control_parameters
1486
1487    CHARACTER(LEN=16) ::  spec_name
1488    CHARACTER(LEN=*)  ::  variable    !<
1489
1490    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
1491    INTEGER(iwp) ::  lsp
1492    INTEGER(iwp) ::  i               !< grid index along x-direction
1493    INTEGER(iwp) ::  j               !< grid index along y-direction
1494    INTEGER(iwp) ::  k               !< grid index along z-direction
1495    INTEGER(iwp) ::  im              !< loop index for masked variables
1496    INTEGER(iwp) ::  jm              !< loop index for masked variables
1497    INTEGER(iwp) ::  kk              !< masked output index along z-direction
1498    INTEGER(iwp) ::  mid             !< masked output running index
1499    INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
1500   
1501    LOGICAL ::  found
1502   
1503    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
1504              local_pf   !<
1505
1506    REAL(wp), PARAMETER ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
1507
1508!
1509!-- local variables.
1510
1511    spec_name = TRIM( variable(4:) )
1512    found = .FALSE.
1513
1514    DO  lsp=1,nspec
1515       IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name) )  THEN             
1516!
1517!-- todo: remove or replace by "CALL message" mechanism (kanani)
1518!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM( variable )  // &
1519!                                                        TRIM( chem_species(lsp)%name )       
1520          IF (av == 0)  THEN
1521             IF ( .NOT. mask_surface(mid) )  THEN
1522
1523                DO  i = 1, mask_size_l(mid,1)
1524                   DO  j = 1, mask_size_l(mid,2) 
1525                      DO  k = 1, mask_size(mid,3) 
1526                          local_pf(i,j,k) = chem_species(lsp)%conc(  &
1527                                               mask_k(mid,k),        &
1528                                               mask_j(mid,j),        &
1529                                               mask_i(mid,i)      )
1530                      ENDDO
1531                   ENDDO
1532                ENDDO
1533
1534             ELSE
1535!             
1536!--             Terrain-following masked output
1537                DO  i = 1, mask_size_l(mid,1)
1538                   DO  j = 1, mask_size_l(mid,2)
1539!--                   Get k index of the highest terraing surface
1540                      im = mask_i(mid,i)
1541                      jm = mask_j(mid,j)
1542                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), DIM = 1 ) - 1
1543                      DO  k = 1, mask_size_l(mid,3)
1544                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
1545!--                      Set value if not in building
1546                         IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
1547                            local_pf(i,j,k) = fill_value
1548                         ELSE
1549                            local_pf(i,j,k) = chem_species(lsp)%conc(kk,jm,im)
1550                         ENDIF
1551                      ENDDO
1552                   ENDDO
1553                ENDDO
1554
1555             ENDIF
1556          ELSE
1557             IF ( .NOT. mask_surface(mid) )  THEN
1558
1559                DO  i = 1, mask_size_l(mid,1)
1560                   DO  j = 1, mask_size_l(mid,2)
1561                      DO  k =  1, mask_size_l(mid,3)
1562                          local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
1563                                               mask_k(mid,k),           &
1564                                               mask_j(mid,j),           &
1565                                               mask_i(mid,i)         )
1566                      ENDDO
1567                   ENDDO
1568                ENDDO
1569
1570             ELSE
1571!             
1572!--             Terrain-following masked output
1573                DO  i = 1, mask_size_l(mid,1)
1574                   DO  j = 1, mask_size_l(mid,2)
1575!--                   Get k index of the highest terraing surface
1576                      im = mask_i(mid,i)
1577                      jm = mask_j(mid,j)
1578                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), DIM = 1 ) - 1
1579                      DO  k = 1, mask_size_l(mid,3)
1580                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
1581!--                      Set value if not in building
1582                         IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
1583                            local_pf(i,j,k) = fill_value
1584                         ELSE
1585                            local_pf(i,j,k) = chem_species(lsp)%conc_av(kk,jm,im)
1586                         ENDIF
1587                      ENDDO
1588                   ENDDO
1589                ENDDO
1590
1591             ENDIF
1592
1593          ENDIF
1594          found = .TRUE.
1595          EXIT
1596       ENDIF
1597    ENDDO
1598
1599    RETURN
1600
1601 END SUBROUTINE chem_data_output_mask     
1602
1603
1604!------------------------------------------------------------------------------!
1605! Description:
1606! ------------
1607!> Subroutine defining appropriate grid for netcdf variables.
1608!> It is called out from subroutine netcdf.
1609!------------------------------------------------------------------------------!
1610 SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1611
1612
1613    CHARACTER (LEN=*), INTENT(IN)  ::  var          !<
1614    LOGICAL, INTENT(OUT)           ::  found        !<
1615    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x       !<
1616    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y       !<
1617    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z       !<
1618
1619    found  = .TRUE.
1620
1621    IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )  THEN                   !< always the same grid for chemistry variables
1622       grid_x = 'x'
1623       grid_y = 'y'
1624       grid_z = 'zu'                             
1625    ELSE
1626       found  = .FALSE.
1627       grid_x = 'none'
1628       grid_y = 'none'
1629       grid_z = 'none'
1630    ENDIF
1631
1632
1633 END SUBROUTINE chem_define_netcdf_grid
1634
1635
1636!------------------------------------------------------------------------------!
1637! Description:
1638! ------------
1639!> Subroutine defining header output for chemistry model
1640!------------------------------------------------------------------------------!
1641 SUBROUTINE chem_header( io )
1642
1643
1644    INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1645    INTEGER(iwp)  :: lsp                       !< running index for chem spcs
1646    INTEGER(iwp)  :: cs_fixed
1647    CHARACTER (LEN=80)  :: docsflux_chr
1648    CHARACTER (LEN=80)  :: docsinit_chr
1649!
1650! Get name of chemical mechanism from chem_gasphase_mod
1651    CALL get_mechanism_name
1652!
1653!-- Write chemistry model  header
1654    WRITE( io, 1 )
1655!
1656!-- Gasphase reaction status
1657    IF ( chem_gasphase_on )  THEN
1658       WRITE( io, 2 )
1659    ELSE
1660       WRITE( io, 3 )
1661    ENDIF
1662!
1663!-- Chemistry time-step
1664    WRITE ( io, 4 ) cs_time_step
1665!
1666!-- Emission mode info
1667!-- At the moment the evaluation is done with both emiss_lod and mode_emis
1668!-- but once salsa has been migrated to emiss_lod the .OR. mode_emis
1669!-- conditions can be removed (ecc 20190513)
1670    IF     ( (emiss_lod == 1) .OR. (mode_emis == 'DEFAULT') )        THEN
1671        WRITE ( io, 5 )
1672    ELSEIF ( (emiss_lod == 0) .OR. (mode_emis == 'PARAMETERIZED') )  THEN
1673        WRITE ( io, 6 )
1674    ELSEIF ( (emiss_lod == 2) .OR. (mode_emis == 'PRE-PROCESSED') )  THEN
1675        WRITE ( io, 7 )
1676    ENDIF
1677!
1678!-- Photolysis scheme info
1679    IF ( photolysis_scheme == "simple" )  THEN
1680       WRITE( io, 8 ) 
1681    ELSEIF (photolysis_scheme == "constant" )  THEN
1682       WRITE( io, 9 )
1683    ENDIF
1684!
1685!-- Emission flux info
1686    lsp = 1
1687    docsflux_chr ='Chemical species for surface emission flux: ' 
1688    DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
1689       docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' 
1690       IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
1691          WRITE ( io, 10 )  docsflux_chr
1692          docsflux_chr = '       '
1693       ENDIF
1694       lsp = lsp + 1
1695    ENDDO
1696
1697    IF ( docsflux_chr /= '' )  THEN
1698       WRITE ( io, 10 )  docsflux_chr
1699    ENDIF
1700!
1701!-- initializatoin of Surface and profile chemical species
1702
1703    lsp = 1
1704    docsinit_chr ='Chemical species for initial surface and profile emissions: ' 
1705    DO WHILE ( cs_name(lsp) /= 'novalue' )
1706       docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ',' 
1707       IF ( LEN_TRIM( docsinit_chr ) >= 75 )  THEN
1708        WRITE ( io, 11 )  docsinit_chr
1709        docsinit_chr = '       '
1710       ENDIF
1711       lsp = lsp + 1
1712    ENDDO
1713
1714    IF ( docsinit_chr /= '' )  THEN
1715       WRITE ( io, 11 )  docsinit_chr
1716    ENDIF
1717
1718    IF ( nesting_chem )  WRITE( io, 12 )  nesting_chem
1719    IF ( nesting_offline_chem )  WRITE( io, 13 )  nesting_offline_chem
1720!
1721!-- number of variable and fix chemical species and number of reactions
1722    cs_fixed = nspec - nvar
1723
1724    WRITE ( io, * ) '   --> Chemical Mechanism        : ', cs_mech
1725    WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
1726    WRITE ( io, * ) '   --> Chemical species, fixed   : ', cs_fixed
1727    WRITE ( io, * ) '   --> Total number of reactions : ', nreact
1728
1729
17301   FORMAT (//' Chemistry model information:'/                                  &
1731           ' ----------------------------'/)
17322   FORMAT ('    --> Chemical reactions are turned on')
17333   FORMAT ('    --> Chemical reactions are turned off')
17344   FORMAT ('    --> Time-step for chemical species: ',F6.2, ' s')
17355   FORMAT ('    --> Emission mode = DEFAULT ')
17366   FORMAT ('    --> Emission mode = PARAMETERIZED ')
17377   FORMAT ('    --> Emission mode = PRE-PROCESSED ')
17388   FORMAT ('    --> Photolysis scheme used =  simple ')
17399   FORMAT ('    --> Photolysis scheme used =  constant ')
174010  FORMAT (/'    ',A) 
174111  FORMAT (/'    ',A) 
174212  FORMAT (/'   Nesting for chemistry variables: ', L1 )
174313  FORMAT (/'   Offline nesting for chemistry variables: ', L1 )
1744!
1745!
1746 END SUBROUTINE chem_header
1747
1748
1749!------------------------------------------------------------------------------!
1750! Description:
1751! ------------
1752!> Subroutine initializating chemistry_model_mod specific arrays
1753!------------------------------------------------------------------------------!
1754 SUBROUTINE chem_init_arrays
1755!
1756!-- Please use this place to allocate required arrays
1757
1758 END SUBROUTINE chem_init_arrays
1759
1760
1761!------------------------------------------------------------------------------!
1762! Description:
1763! ------------
1764!> Subroutine initializating chemistry_model_mod
1765!------------------------------------------------------------------------------!
1766 SUBROUTINE chem_init
1767
1768    USE chem_emissions_mod,                                                    &
1769        ONLY:  chem_emissions_init
1770       
1771    USE netcdf_data_input_mod,                                                 &
1772        ONLY:  init_3d
1773
1774
1775    INTEGER(iwp) ::  i !< running index x dimension
1776    INTEGER(iwp) ::  j !< running index y dimension
1777    INTEGER(iwp) ::  n !< running index for chemical species
1778
1779
1780    IF ( debug_output )  CALL debug_message( 'chem_init', 'start' )
1781!
1782!-- Next statement is to avoid compiler warning about unused variables
1783    IF ( ( ilu_arable + ilu_coniferous_forest + ilu_deciduous_forest + ilu_mediterrean_scrub + &
1784           ilu_permanent_crops + ilu_savanna + ilu_semi_natural_veg + ilu_tropical_forest +    &
1785           ilu_urban ) == 0 )  CONTINUE
1786         
1787    IF ( emissions_anthropogenic )  CALL chem_emissions_init
1788!
1789!-- Chemistry variables will be initialized if availabe from dynamic
1790!-- input file. Note, it is possible to initialize only part of the chemistry
1791!-- variables from dynamic input.
1792    IF ( INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
1793       DO  n = 1, nspec
1794          IF ( init_3d%from_file_chem(n) )  THEN
1795             DO  i = nxlg, nxrg
1796                DO  j = nysg, nyng
1797                   chem_species(n)%conc(:,j,i) = init_3d%chem_init(:,n)
1798                ENDDO
1799             ENDDO
1800          ENDIF
1801       ENDDO
1802    ENDIF
1803
1804    IF ( debug_output )  CALL debug_message( 'chem_init', 'end' )
1805
1806 END SUBROUTINE chem_init
1807
1808
1809!------------------------------------------------------------------------------!
1810! Description:
1811! ------------
1812!> Subroutine initializating chemistry_model_mod
1813!> internal workaround for chem_species dependency in chem_check_parameters
1814!------------------------------------------------------------------------------!
1815 SUBROUTINE chem_init_internal
1816 
1817    USE pegrid
1818
1819    USE netcdf_data_input_mod,                                                 &
1820        ONLY:  chem_emis, chem_emis_att, input_pids_dynamic, init_3d,          &
1821               netcdf_data_input_chemistry_data
1822
1823!
1824!-- Local variables
1825    INTEGER(iwp) ::  i                 !< running index for for horiz numerical grid points
1826    INTEGER(iwp) ::  j                 !< running index for for horiz numerical grid points
1827    INTEGER(iwp) ::  lsp               !< running index for chem spcs
1828    INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
1829
1830    IF ( emissions_anthropogenic )  THEN
1831       CALL netcdf_data_input_chemistry_data( chem_emis_att, chem_emis )
1832    ENDIF
1833!
1834!-- Allocate memory for chemical species
1835    ALLOCATE( chem_species(nspec) )
1836    ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1837    ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1838    ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1839    ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) 
1840    ALLOCATE( phot_frequen(nphot) ) 
1841    ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
1842    ALLOCATE( bc_cs_t_val(nspec) )
1843!
1844!-- Initialize arrays
1845    spec_conc_1 (:,:,:,:) = 0.0_wp
1846    spec_conc_2 (:,:,:,:) = 0.0_wp
1847    spec_conc_3 (:,:,:,:) = 0.0_wp
1848    spec_conc_av(:,:,:,:) = 0.0_wp
1849
1850
1851    DO  lsp = 1, nspec
1852       chem_species(lsp)%name    = spc_names(lsp)
1853
1854       chem_species(lsp)%conc   (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_1 (:,:,:,lsp)
1855       chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_2 (:,:,:,lsp)
1856       chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
1857       chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_av(:,:,:,lsp)     
1858
1859       ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))                   
1860       chem_species(lsp)%cssws_av    = 0.0_wp
1861!
1862!--    The following block can be useful when emission module is not applied. &
1863!--    if emission module is applied the following block will be overwritten.
1864       ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))   
1865       ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))   
1866       ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) 
1867       ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))   
1868       chem_species(lsp)%flux_s_cs = 0.0_wp                                     
1869       chem_species(lsp)%flux_l_cs = 0.0_wp                                     
1870       chem_species(lsp)%diss_s_cs = 0.0_wp                                     
1871       chem_species(lsp)%diss_l_cs = 0.0_wp                                     
1872!
1873!--   Allocate memory for initial concentration profiles
1874!--   (concentration values come from namelist)
1875!--   (@todo (FK): Because of this, chem_init is called in palm before
1876!--               check_parameters, since conc_pr_init is used there.
1877!--               We have to find another solution since chem_init should
1878!--               eventually be called from init_3d_model!!)
1879       ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
1880       chem_species(lsp)%conc_pr_init(:) = 0.0_wp
1881
1882    ENDDO
1883!
1884!-- Set control flags for decycling only at lateral boundary cores, within the
1885!-- inner cores the decycle flag is set to .False.. Even though it does not
1886!-- affect the setting of chemistry boundary conditions, this flag is used to
1887!-- set advection control flags appropriately.
1888    decycle_chem_lr = MERGE( decycle_chem_lr, .FALSE.,                         &
1889                             nxl == 0  .OR.  nxr == nx )
1890    decycle_chem_ns = MERGE( decycle_chem_ns, .FALSE.,                         &
1891                             nys == 0  .OR.  nyn == ny )
1892!
1893!-- For some passive scalars decycling may be enabled. This case, the lateral
1894!-- boundary conditions are non-cyclic for these scalars (chemical species
1895!-- and aerosols), while the other scalars may have
1896!-- cyclic boundary conditions. However, large gradients near the boundaries
1897!-- may produce stationary numerical oscillations near the lateral boundaries
1898!-- when a higher-order scheme is applied near these boundaries.
1899!-- To get rid-off this, set-up additional flags that control the order of the
1900!-- scalar advection scheme near the lateral boundaries for passive scalars
1901!-- with decycling.
1902    IF ( scalar_advec == 'ws-scheme' )  THEN
1903       ALLOCATE( cs_advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1904!
1905!--    In case of decyling, set Neumann boundary conditions for wall_flags_total_0
1906!--    bit 31 instead of cyclic boundary conditions.
1907!--    Bit 31 is used to identify extended degradation zones (please see
1908!--    following comment).
1909!--    Note, since several also other modules like Salsa or other future
1910!--    one may access this bit but may have other boundary conditions, the
1911!--    original value of wall_flags_total_0 bit 31 must not be modified. Hence,
1912!--    store the boundary conditions directly on cs_advc_flags_s.
1913!--    cs_advc_flags_s will be later overwritten in ws_init_flags_scalar and
1914!--    bit 31 won't be used to control the numerical order.
1915!--    Initialize with flag 31 only.
1916       cs_advc_flags_s = 0
1917       cs_advc_flags_s = MERGE( IBSET( cs_advc_flags_s, 31 ), 0,               &
1918                                BTEST( wall_flags_total_0, 31 ) )
1919
1920       IF ( decycle_chem_ns )  THEN
1921          IF ( nys == 0  )  THEN
1922             DO  i = 1, nbgp     
1923                cs_advc_flags_s(:,nys-i,:) = MERGE(                            &
1924                                        IBSET( cs_advc_flags_s(:,nys,:), 31 ), &
1925                                        IBCLR( cs_advc_flags_s(:,nys,:), 31 ), &
1926                                        BTEST( cs_advc_flags_s(:,nys,:), 31 )  &
1927                                                  )
1928             ENDDO
1929          ENDIF
1930          IF ( nyn == ny )  THEN
1931             DO  i = 1, nbgp 
1932                cs_advc_flags_s(:,nyn+i,:) = MERGE(                            &
1933                                        IBSET( cs_advc_flags_s(:,nyn,:), 31 ), &
1934                                        IBCLR( cs_advc_flags_s(:,nyn,:), 31 ), &
1935                                        BTEST( cs_advc_flags_s(:,nyn,:), 31 )  &
1936                                                  )
1937             ENDDO
1938          ENDIF
1939       ENDIF
1940       IF ( decycle_chem_lr )  THEN
1941          IF ( nxl == 0  )  THEN
1942             DO  i = 1, nbgp   
1943                cs_advc_flags_s(:,:,nxl-i) = MERGE(                            &
1944                                        IBSET( cs_advc_flags_s(:,:,nxl), 31 ), &
1945                                        IBCLR( cs_advc_flags_s(:,:,nxl), 31 ), &
1946                                        BTEST( cs_advc_flags_s(:,:,nxl), 31 )  &
1947                                                  )
1948             ENDDO
1949          ENDIF
1950          IF ( nxr == nx )  THEN
1951             DO  i = 1, nbgp   
1952                cs_advc_flags_s(:,:,nxr+i) = MERGE(                            &
1953                                        IBSET( cs_advc_flags_s(:,:,nxr), 31 ), &
1954                                        IBCLR( cs_advc_flags_s(:,:,nxr), 31 ), &
1955                                        BTEST( cs_advc_flags_s(:,:,nxr), 31 )  &
1956                                                  )
1957             ENDDO
1958          ENDIF     
1959       ENDIF
1960!
1961!--    To initialize advection flags appropriately, pass the boundary flags.
1962!--    The last argument indicates that a passive scalar is treated, where
1963!--    the horizontal advection terms are degraded already 2 grid points before
1964!--    the lateral boundary to avoid stationary oscillations at large-gradients.
1965!--    Also, extended degradation zones are applied, where horizontal advection of
1966!--    passive scalars is discretized by first-order scheme at all grid points
1967!--    that in the vicinity of buildings (<= 3 grid points). Even though no
1968!--    building is within the numerical stencil, first-order scheme is used.
1969!--    At fourth and fifth grid point the order of the horizontal advection scheme
1970!--    is successively upgraded.
1971!--    These extended degradation zones are used to avoid stationary numerical
1972!--    oscillations, which are responsible for high concentration maxima that may
1973!--    appear under shear-free stable conditions.
1974       CALL ws_init_flags_scalar(                                              &
1975                  bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_chem_lr, &
1976                  bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_chem_ns, &
1977                  bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_chem_lr, &
1978                  bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_chem_ns, &
1979                  cs_advc_flags_s, .TRUE. )
1980    ENDIF
1981!
1982!-- Initial concentration of profiles is prescribed by parameters cs_profile
1983!-- and cs_heights in the namelist &chemistry_parameters
1984
1985    CALL chem_init_profiles
1986!   
1987!-- In case there is dynamic input file, create a list of names for chemistry
1988!-- initial input files. Also, initialize array that indicates whether the
1989!-- respective variable is on file or not.
1990    IF ( input_pids_dynamic )  THEN   
1991       ALLOCATE( init_3d%var_names_chem(1:nspec) )
1992       ALLOCATE( init_3d%from_file_chem(1:nspec) )
1993       init_3d%from_file_chem(:) = .FALSE.
1994       
1995       DO  lsp = 1, nspec
1996          init_3d%var_names_chem(lsp) = init_3d%init_char // TRIM( chem_species(lsp)%name )
1997       ENDDO
1998    ENDIF
1999!
2000!-- Initialize model variables
2001    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
2002         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
2003!
2004!--    First model run of a possible job queue.
2005!--    Initial profiles of the variables must be computed.
2006       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
2007!
2008!--       Transfer initial profiles to the arrays of the 3D model
2009          DO  lsp = 1, nspec
2010             DO  i = nxlg, nxrg
2011                DO  j = nysg, nyng
2012                   DO lpr_lev = 1, nz + 1
2013                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
2014                   ENDDO
2015                ENDDO
2016             ENDDO   
2017          ENDDO
2018
2019       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
2020       THEN
2021
2022          DO  lsp = 1, nspec
2023             DO  i = nxlg, nxrg
2024                DO  j = nysg, nyng
2025                   chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init   
2026                ENDDO
2027             ENDDO
2028          ENDDO
2029
2030       ENDIF
2031!
2032!--    If required, change the surface chem spcs at the start of the 3D run
2033       IF ( cs_surface_initial_change(1) /= 0.0_wp )  THEN           
2034          DO  lsp = 1, nspec
2035             chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
2036                                               cs_surface_initial_change(lsp)
2037          ENDDO
2038       ENDIF
2039
2040    ENDIF
2041!
2042!-- Initiale old and new time levels. Note, this has to be done also in restart runs
2043    DO  lsp = 1, nvar
2044       chem_species(lsp)%tconc_m = 0.0_wp                     
2045       chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
2046    ENDDO
2047
2048    DO  lsp = 1, nphot
2049       phot_frequen(lsp)%name = phot_names(lsp)
2050!
2051!-- todo: remove or replace by "CALL message" mechanism (kanani)
2052!--       IF( myid == 0 )  THEN
2053!--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
2054!--       ENDIF
2055       phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
2056    ENDDO
2057
2058!    CALL photolysis_init   ! probably also required for restart
2059
2060    RETURN
2061
2062 END SUBROUTINE chem_init_internal
2063
2064
2065!------------------------------------------------------------------------------!
2066! Description:
2067! ------------
2068!> Subroutine defining initial vertical profiles of chemical species (given by
2069!> namelist parameters chem_profiles and chem_heights)  --> which should work
2070!> analogue to parameters u_profile, v_profile and uv_heights)
2071!------------------------------------------------------------------------------!
2072 SUBROUTINE chem_init_profiles             
2073
2074    USE chem_modules
2075
2076!
2077!-- Local variables
2078    INTEGER ::  lsp        !< running index for number of species in derived data type species_def
2079    INTEGER ::  lsp_usr     !< running index for number of species (user defined)  in cs_names, cs_profiles etc
2080    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
2081    INTEGER ::  npr_lev    !< the next available profile lev
2082!
2083!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
2084!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
2085!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
2086!-- "cs_heights" are prescribed, their values will!override the constant profile given by
2087!-- "cs_surface".
2088!     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2089    lsp_usr = 1
2090    DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
2091       DO  lsp = 1, nspec                                !
2092!
2093!--       create initial profile (conc_pr_init) for each chemical species
2094          IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
2095             IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
2096!
2097!--            set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
2098                DO lpr_lev = 0, nzt+1
2099                   chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
2100                ENDDO
2101             ELSE
2102                IF ( cs_heights(1,1) /= 0.0_wp )  THEN
2103                   message_string = 'The surface value of cs_heights must be 0.0'
2104                   CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
2105                ENDIF
2106
2107                use_prescribed_profile_data = .TRUE.
2108
2109                npr_lev = 1
2110!                chem_species(lsp)%conc_pr_init(0) = 0.0_wp
2111                DO  lpr_lev = 1, nz+1
2112                   IF ( npr_lev < 100 )  THEN
2113                      DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
2114                         npr_lev = npr_lev + 1
2115                         IF ( npr_lev == 100 )  THEN
2116                            message_string = 'number of chem spcs exceeding the limit'
2117                            CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
2118                            EXIT
2119                         ENDIF
2120                      ENDDO
2121                   ENDIF
2122                   IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
2123                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
2124                           ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
2125                           ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
2126                           ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
2127                   ELSE
2128                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
2129                   ENDIF
2130                ENDDO
2131             ENDIF
2132!
2133!--       If a profile is prescribed explicity using cs_profiles and cs_heights, then 
2134!--       chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
2135!--       on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
2136          ENDIF
2137
2138       ENDDO
2139
2140       lsp_usr = lsp_usr + 1
2141    ENDDO
2142!     ENDIF
2143
2144 END SUBROUTINE chem_init_profiles
2145
2146 
2147!------------------------------------------------------------------------------!
2148! Description:
2149! ------------
2150!> Subroutine to integrate chemical species in the given chemical mechanism
2151!------------------------------------------------------------------------------!
2152 SUBROUTINE chem_integrate_ij( i, j )
2153
2154    USE statistics,                                                          &
2155        ONLY:  weight_pres
2156
2157    USE control_parameters,                                                  &
2158        ONLY:  dt_3d, intermediate_timestep_count, time_since_reference_point
2159
2160
2161    INTEGER,INTENT(IN)       :: i
2162    INTEGER,INTENT(IN)       :: j
2163!
2164!--   local variables
2165    INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
2166    INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
2167    INTEGER, DIMENSION(20)    :: istatus
2168    REAL(kind=wp), DIMENSION(nzb+1:nzt,nspec)                :: tmp_conc
2169    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_temp
2170    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_qvap
2171    REAL(kind=wp), DIMENSION(nzb+1:nzt,nphot)                :: tmp_phot
2172    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact
2173    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
2174                                                                              !<    molecules cm^{-3} and ppm
2175
2176    INTEGER,DIMENSION(nzb+1:nzt)                            :: nacc          !< Number of accepted steps
2177    INTEGER,DIMENSION(nzb+1:nzt)                            :: nrej          !< Number of rejected steps
2178
2179    REAL(wp)                         ::  conv                                !< conversion factor
2180    REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp                 !< Conversion factor ppm to fraction
2181    REAL(wp), PARAMETER              ::  fr2ppm  = 1.0e6_wp                  !< Conversion factor fraction to ppm
2182!    REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp                  !< Mole mass of dry air
2183!    REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp               !< Mole mass of water vapor
2184    REAL(wp), PARAMETER              ::  t_std   = 273.15_wp                 !< standard pressure (Pa)
2185    REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp               !< standard pressure (Pa)
2186    REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp               !< Mole volume (22.414 l) in cm^3
2187    REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp               !< Avogadro number (molecules/mol)
2188
2189    REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
2190
2191    REAL(kind=wp)  :: dt_chem                                             
2192!
2193!-- Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
2194    IF (chem_gasphase_on)  THEN
2195       nacc = 0
2196       nrej = 0
2197
2198       tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
2199!
2200!--    convert ppm to molecules/cm**3
2201!--    tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp *
2202!--               hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 
2203       conv = ppm2fr * xna / vmolcm
2204       tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
2205       tmp_fact_i = 1.0_wp/tmp_fact
2206
2207       IF ( humidity )  THEN
2208          IF ( bulk_cloud_model )  THEN
2209             tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) *  &
2210                             xm_air/xm_h2o * fr2ppm * tmp_fact(:)
2211          ELSE
2212             tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * fr2ppm * tmp_fact(:)
2213          ENDIF
2214       ELSE
2215          tmp_qvap(:) = 0.01 * xm_air/xm_h2o * fr2ppm * tmp_fact(:)          !< Constant value for q if water vapor is not computed
2216       ENDIF
2217
2218       DO  lsp = 1,nspec
2219          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) 
2220       ENDDO
2221
2222       DO lph = 1,nphot
2223          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
2224       ENDDO
2225!
2226!--    Compute length of time step
2227       IF ( call_chem_at_all_substeps )  THEN
2228          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
2229       ELSE
2230          dt_chem = dt_3d
2231       ENDIF
2232
2233       cs_time_step = dt_chem
2234
2235       IF(maxval(rcntrl) > 0.0)  THEN    ! Only if rcntrl is set
2236          IF( time_since_reference_point <= 2*dt_3d)  THEN
2237             rcntrl_local = 0
2238          ELSE
2239             rcntrl_local = rcntrl
2240          ENDIF
2241       ELSE
2242          rcntrl_local = 0
2243       END IF
2244
2245       CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
2246            icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
2247
2248       DO  lsp = 1,nspec
2249          chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
2250       ENDDO
2251
2252
2253    ENDIF
2254
2255    RETURN
2256 END SUBROUTINE chem_integrate_ij
2257
2258
2259!------------------------------------------------------------------------------!
2260! Description:
2261! ------------
2262!> Subroutine defining parin for &chemistry_parameters for chemistry model
2263!------------------------------------------------------------------------------!
2264 SUBROUTINE chem_parin
2265
2266    USE chem_modules
2267    USE control_parameters
2268
2269    USE pegrid
2270    USE statistics
2271
2272
2273    CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current line of the parameter file
2274
2275    REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0 automatic stepping
2276    INTEGER(iwp) ::  i                                 !<
2277    INTEGER(iwp) ::  max_pr_cs_tmp                     !<
2278
2279
2280    NAMELIST /chemistry_parameters/  bc_cs_b,                          &
2281         bc_cs_t,                          &
2282         call_chem_at_all_substeps,        &
2283         chem_debug0,                      &
2284         chem_debug1,                      &
2285         chem_debug2,                      &
2286         chem_gasphase_on,                 &
2287         chem_mechanism,                   &         
2288         cs_heights,                       &
2289         cs_name,                          &
2290         cs_profile,                       &
2291         cs_surface,                       &
2292         cs_surface_initial_change,        &
2293         cs_vertical_gradient_level,       &
2294         daytype_mdh,                      &
2295         decycle_chem_lr,                  &
2296         decycle_chem_ns,                  &           
2297         decycle_method,                   &
2298         deposition_dry,                   &
2299         emissions_anthropogenic,          &
2300         emiss_lod,                        &
2301         emiss_factor_main,                &
2302         emiss_factor_side,                &                     
2303         icntrl,                           &
2304         main_street_id,                   &
2305         max_street_id,                    &
2306         mode_emis,                        &
2307         my_steps,                         &
2308         nesting_chem,                     &
2309         nesting_offline_chem,             &
2310         rcntrl,                           &
2311         side_street_id,                   &
2312         photolysis_scheme,                &
2313         wall_csflux,                      &
2314         cs_vertical_gradient,             &
2315         top_csflux,                       & 
2316         surface_csflux,                   &
2317         surface_csflux_name,              &
2318         time_fac_type
2319!
2320!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
2321!-- so this way we could prescribe a specific flux value for each species
2322    !>  chemistry_parameters for initial profiles
2323    !>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
2324    !>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
2325    !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
2326    !>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
2327    !>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
2328    !>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
2329!
2330!--   Read chem namelist   
2331
2332    CHARACTER(LEN=8)    :: solver_type
2333
2334    icntrl    = 0
2335    rcntrl    = 0.0_wp
2336    my_steps  = 0.0_wp
2337    photolysis_scheme = 'simple'
2338    atol = 1.0_wp
2339    rtol = 0.01_wp
2340!
2341!--   Try to find chemistry package
2342    REWIND ( 11 )
2343    line = ' '
2344    DO   WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
2345       READ ( 11, '(A)', END=20 )  line
2346    ENDDO
2347    BACKSPACE ( 11 )
2348!
2349!--   Read chemistry namelist
2350    READ ( 11, chemistry_parameters, ERR = 10, END = 20 )     
2351!
2352!--   Enable chemistry model
2353    air_chemistry = .TRUE.                   
2354    GOTO 20
2355
2356 10 BACKSPACE( 11 )
2357    READ( 11 , '(A)') line
2358    CALL parin_fail_message( 'chemistry_parameters', line )
2359
2360 20 CONTINUE
2361
2362!
2363!-- synchronize emiss_lod and mod_emis only if emissions_anthropogenic
2364!-- is activated in the namelist.  Otherwise their values are "don't care"
2365    IF ( emissions_anthropogenic )  THEN
2366!
2367!--    check for emission mode for chem species
2368
2369       IF ( emiss_lod < 0 )  THEN   !- if LOD not defined in namelist
2370          IF ( ( mode_emis /= 'PARAMETERIZED'  )    .AND.      &
2371               ( mode_emis /= 'DEFAULT'        )    .AND.      &
2372               ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
2373             message_string = 'Incorrect mode_emiss  option select. Please check spelling'
2374             CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2375          ENDIF
2376       ELSE
2377          IF ( ( emiss_lod /= 0 )    .AND.         &
2378               ( emiss_lod /= 1 )    .AND.         &
2379               ( emiss_lod /= 2 ) )  THEN
2380             message_string = 'Invalid value for emiss_lod (0, 1, or 2)'
2381             CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2382          ENDIF
2383       ENDIF
2384
2385!
2386! for reference (ecc)
2387!    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
2388!       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
2389!       CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
2390!    ENDIF
2391
2392!
2393!-- conflict resolution for emiss_lod and mode_emis
2394!-- 1) if emiss_lod is defined, have mode_emis assume same setting as emiss_lod
2395!-- 2) if emiss_lod it not defined, have emiss_lod assuem same setting as mode_emis
2396!-- this check is in place to retain backward compatibility with salsa until the
2397!-- code is migrated completed to emiss_lod
2398!-- note that
2399
2400       IF  ( emiss_lod >= 0 ) THEN
2401
2402          SELECT CASE  ( emiss_lod )
2403             CASE (0)  !- parameterized mode
2404                mode_emis = 'PARAMETERIZED'
2405             CASE (1)  !- default mode
2406                mode_emis = 'DEFAULT'
2407             CASE (2)  !- preprocessed mode
2408                mode_emis = 'PRE-PROCESSED'
2409          END SELECT
2410       
2411          message_string = 'Synchronizing mode_emis to defined emiss_lod'               //  &
2412                           CHAR(10)  //  '          '                                   //  &
2413                           'NOTE - mode_emis will be depreciated in future releases'    //  &
2414                           CHAR(10)  //  '          '                                   //  &
2415                           'please use emiss_lod to define emission mode'
2416 
2417          CALL message ( 'parin_chem', 'CM0463', 0, 0, 0, 6, 0 )
2418
2419       ELSE ! if emiss_lod is not set
2420
2421          SELECT CASE ( mode_emis )
2422             CASE ('PARAMETERIZED')
2423                emiss_lod = 0
2424             CASE ('DEFAULT')
2425                emiss_lod = 1
2426             CASE ('PRE-PROCESSED')
2427                emiss_lod = 2
2428          END SELECT
2429
2430          message_string = 'emiss_lod undefined.  Using existing mod_emiss setting'     //  &
2431                           CHAR(10)  //  '          '                                   //  &
2432                           'NOTE - mode_emis will be depreciated in future releases'    //  &
2433                           CHAR(10)  //  '          '                                   //  &
2434                           '       please use emiss_lod to define emission mode'
2435
2436          CALL message ( 'parin_chem', 'CM0464', 0, 0, 0, 6, 0 )
2437       ENDIF
2438
2439    ENDIF  ! if emissions_anthropengic
2440
2441    t_steps = my_steps         
2442!
2443!--    Determine the number of user-defined profiles and append them to the
2444!--    standard data output (data_output_pr)
2445    max_pr_cs_tmp = 0 
2446    i = 1
2447
2448    DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
2449       IF ( TRIM( data_output_pr(i)(1:3) ) == 'kc_' )  THEN
2450          max_pr_cs_tmp = max_pr_cs_tmp+1
2451       ENDIF
2452       i = i +1
2453    ENDDO
2454
2455    IF ( max_pr_cs_tmp > 0 )  THEN
2456       cs_pr_namelist_found = .TRUE.
2457       max_pr_cs = max_pr_cs_tmp
2458    ENDIF
2459
2460    !      Set Solver Type
2461    IF(icntrl(3) == 0)  THEN
2462       solver_type = 'rodas3'           !Default
2463    ELSE IF(icntrl(3) == 1)  THEN
2464       solver_type = 'ros2'
2465    ELSE IF(icntrl(3) == 2)  THEN
2466       solver_type = 'ros3'
2467    ELSE IF(icntrl(3) == 3)  THEN
2468       solver_type = 'ro4'
2469    ELSE IF(icntrl(3) == 4)  THEN
2470       solver_type = 'rodas3'
2471    ELSE IF(icntrl(3) == 5)  THEN
2472       solver_type = 'rodas4'
2473    ELSE IF(icntrl(3) == 6)  THEN
2474       solver_type = 'Rang3'
2475    ELSE
2476       message_string = 'illegal solver type'
2477       CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
2478    END IF
2479
2480!
2481!--   todo: remove or replace by "CALL message" mechanism (kanani)
2482!       write(text,*) 'gas_phase chemistry: solver_type = ',TRIM( solver_type )
2483!kk    Has to be changed to right calling sequence
2484!        IF(myid == 0)  THEN
2485!           write(9,*) ' '
2486!           write(9,*) 'kpp setup '
2487!           write(9,*) ' '
2488!           write(9,*) '    gas_phase chemistry: solver_type = ',TRIM( solver_type )
2489!           write(9,*) ' '
2490!           write(9,*) '    Hstart  = ',rcntrl(3)
2491!           write(9,*) '    FacMin  = ',rcntrl(4)
2492!           write(9,*) '    FacMax  = ',rcntrl(5)
2493!           write(9,*) ' '
2494!           IF(vl_dim > 1)  THEN
2495!              write(9,*) '    Vector mode                   vektor length = ',vl_dim
2496!           ELSE
2497!              write(9,*) '    Scalar mode'
2498!           ENDIF
2499!           write(9,*) ' '
2500!        END IF
2501
2502    RETURN
2503
2504 END SUBROUTINE chem_parin
2505
2506
2507!------------------------------------------------------------------------------!
2508! Description:
2509! ------------
2510!> Call for all grid points
2511!------------------------------------------------------------------------------!
2512    SUBROUTINE chem_actions( location )
2513
2514
2515    CHARACTER (LEN=*), INTENT(IN) ::  location !< call location string
2516
2517    SELECT CASE ( location )
2518
2519       CASE ( 'before_prognostic_equations' )
2520!
2521!--       Chemical reactions and deposition
2522          IF ( chem_gasphase_on ) THEN
2523!
2524!--          If required, calculate photolysis frequencies -
2525!--          UNFINISHED: Why not before the intermediate timestep loop?
2526             IF ( intermediate_timestep_count ==  1 )  THEN
2527                CALL photolysis_control
2528             ENDIF
2529
2530          ENDIF
2531
2532       CASE DEFAULT
2533          CONTINUE
2534
2535    END SELECT
2536
2537    END SUBROUTINE chem_actions
2538
2539
2540!------------------------------------------------------------------------------!
2541! Description:
2542! ------------
2543!> Call for grid points i,j
2544!------------------------------------------------------------------------------!
2545
2546    SUBROUTINE chem_actions_ij( i, j, location )
2547
2548
2549    INTEGER(iwp),      INTENT(IN) ::  i         !< grid index in x-direction
2550    INTEGER(iwp),      INTENT(IN) ::  j         !< grid index in y-direction
2551    CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
2552    INTEGER(iwp)  ::  dummy  !< call location string
2553
2554    IF ( air_chemistry    )   dummy = i + j
2555
2556    SELECT CASE ( location )
2557
2558       CASE DEFAULT
2559          CONTINUE
2560
2561    END SELECT
2562
2563
2564    END SUBROUTINE chem_actions_ij
2565
2566
2567!------------------------------------------------------------------------------!
2568! Description:
2569! ------------
2570!> Call for all grid points
2571!------------------------------------------------------------------------------!
2572    SUBROUTINE chem_non_advective_processes()
2573
2574
2575      INTEGER(iwp) ::  i  !<
2576      INTEGER(iwp) ::  j  !<
2577
2578      !
2579      !-- Calculation of chemical reactions and deposition.
2580
2581
2582      IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
2583
2584         IF ( chem_gasphase_on ) THEN
2585            CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2586            !$OMP PARALLEL PRIVATE (i,j)
2587            !$OMP DO schedule(static,1)
2588            DO  i = nxl, nxr
2589               DO  j = nys, nyn
2590                  CALL chem_integrate( i, j )
2591               ENDDO
2592            ENDDO
2593            !$OMP END PARALLEL
2594            CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2595         ENDIF
2596
2597         IF ( deposition_dry )  THEN
2598            CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2599            DO  i = nxl, nxr
2600               DO  j = nys, nyn
2601                  CALL chem_depo( i, j )
2602               ENDDO
2603            ENDDO
2604            CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2605         ENDIF
2606
2607      ENDIF
2608
2609
2610
2611    END SUBROUTINE chem_non_advective_processes
2612
2613
2614!------------------------------------------------------------------------------!
2615! Description:
2616! ------------
2617!> Call for grid points i,j
2618!------------------------------------------------------------------------------!
2619 SUBROUTINE chem_non_advective_processes_ij( i, j )
2620
2621
2622   INTEGER(iwp), INTENT(IN) ::  i  !< grid index in x-direction
2623   INTEGER(iwp), INTENT(IN) ::  j  !< grid index in y-direction
2624
2625!
2626!-- Calculation of chemical reactions and deposition.
2627
2628
2629   IF ( intermediate_timestep_count == 1 .OR. call_chem_at_all_substeps )  THEN
2630
2631      IF ( chem_gasphase_on ) THEN
2632         CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
2633         CALL chem_integrate( i, j )
2634         CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
2635      ENDIF
2636
2637      IF ( deposition_dry )  THEN
2638         CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
2639         CALL chem_depo( i, j )
2640         CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
2641      ENDIF
2642
2643   ENDIF
2644
2645
2646
2647 END SUBROUTINE chem_non_advective_processes_ij
2648 
2649!------------------------------------------------------------------------------!
2650! Description:
2651! ------------
2652!> routine for exchange horiz of chemical quantities 
2653!------------------------------------------------------------------------------!
2654 SUBROUTINE chem_exchange_horiz_bounds
2655 
2656   INTEGER(iwp) ::  lsp       !<
2657   INTEGER(iwp) ::  lsp_usr   !<
2658 
2659!
2660!--    Loop over chemical species       
2661       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
2662       DO  lsp = 1, nvar
2663          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )   
2664          lsp_usr = 1 
2665          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
2666             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
2667!
2668!--             As chem_exchange_horiz_bounds is called at the beginning
2669!--             of prognostic_equations, boundary conditions are set on
2670!--             %conc.
2671                CALL chem_boundary_conds( chem_species(lsp)%conc,              &
2672                                          chem_species(lsp)%conc_pr_init ) 
2673             
2674             ENDIF
2675             lsp_usr = lsp_usr + 1
2676          ENDDO
2677
2678
2679       ENDDO
2680       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
2681
2682
2683 END SUBROUTINE chem_exchange_horiz_bounds
2684
2685 
2686!------------------------------------------------------------------------------!
2687! Description:
2688! ------------
2689!> Subroutine calculating prognostic equations for chemical species
2690!> (vector-optimized).
2691!> Routine is called separately for each chemical species over a loop from
2692!> prognostic_equations.
2693!------------------------------------------------------------------------------!
2694 SUBROUTINE chem_prognostic_equations()
2695
2696
2697    INTEGER ::  i   !< running index
2698    INTEGER ::  j   !< running index
2699    INTEGER ::  k   !< running index
2700
2701    INTEGER(iwp) ::  ilsp   !<
2702
2703
2704    CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2705
2706    DO  ilsp = 1, nvar
2707!
2708!--    Tendency terms for chemical species
2709       tend = 0.0_wp
2710!
2711!--    Advection terms
2712       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2713          IF ( ws_scheme_sca )  THEN
2714             CALL advec_s_ws( cs_advc_flags_s, chem_species(ilsp)%conc, 'kc',                      &
2715                              bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_chem_lr,         &
2716                              bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_chem_ns,         &
2717                              bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_chem_lr,         &
2718                              bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_chem_ns )
2719          ELSE
2720             CALL advec_s_pw( chem_species(ilsp)%conc )
2721          ENDIF
2722       ELSE
2723          CALL advec_s_up( chem_species(ilsp)%conc )
2724       ENDIF
2725!
2726!--    Diffusion terms  (the last three arguments are zero)
2727       CALL diffusion_s( chem_species(ilsp)%conc,                                                  &
2728            surf_def_h(0)%cssws(ilsp,:),                                                           &
2729            surf_def_h(1)%cssws(ilsp,:),                                                           &
2730            surf_def_h(2)%cssws(ilsp,:),                                                           &
2731            surf_lsm_h%cssws(ilsp,:),                                                              &
2732            surf_usm_h%cssws(ilsp,:),                                                              &
2733            surf_def_v(0)%cssws(ilsp,:),                                                           &
2734            surf_def_v(1)%cssws(ilsp,:),                                                           &
2735            surf_def_v(2)%cssws(ilsp,:),                                                           &
2736            surf_def_v(3)%cssws(ilsp,:),                                                           &
2737            surf_lsm_v(0)%cssws(ilsp,:),                                                           &
2738            surf_lsm_v(1)%cssws(ilsp,:),                                                           &
2739            surf_lsm_v(2)%cssws(ilsp,:),                                                           &
2740            surf_lsm_v(3)%cssws(ilsp,:),                                                           &
2741            surf_usm_v(0)%cssws(ilsp,:),                                                           &
2742            surf_usm_v(1)%cssws(ilsp,:),                                                           &
2743            surf_usm_v(2)%cssws(ilsp,:),                                                           &
2744            surf_usm_v(3)%cssws(ilsp,:) )
2745!
2746!--    Prognostic equation for chemical species
2747       DO  i = nxl, nxr
2748          DO  j = nys, nyn
2749             !following directive is required to vectorize on Intel19
2750             !DIR$ IVDEP
2751             DO  k = nzb+1, nzt
2752                chem_species(ilsp)%conc_p(k,j,i) =   chem_species(ilsp)%conc(k,j,i)                &
2753                     + ( dt_3d  *                                                                  &
2754                     (   tsc(2) * tend(k,j,i)                                                      &
2755                     + tsc(3) * chem_species(ilsp)%tconc_m(k,j,i)                                  &
2756                     )                                                                             &
2757                     - tsc(5) * rdf_sc(k)                                                          &
2758                     * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )     &
2759                     )                                                                             &
2760                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
2761
2762                IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN
2763                   chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)
2764                ENDIF
2765             ENDDO
2766          ENDDO
2767       ENDDO
2768!
2769!--    Calculate tendencies for the next Runge-Kutta step
2770       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2771          IF ( intermediate_timestep_count == 1 )  THEN
2772             DO  i = nxl, nxr
2773                DO  j = nys, nyn
2774                   DO  k = nzb+1, nzt
2775                      chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
2776                   ENDDO
2777                ENDDO
2778             ENDDO
2779          ELSEIF ( intermediate_timestep_count < &
2780               intermediate_timestep_count_max )  THEN
2781             DO  i = nxl, nxr
2782                DO  j = nys, nyn
2783                   DO  k = nzb+1, nzt
2784                      chem_species(ilsp)%tconc_m(k,j,i) = - 9.5625_wp * tend(k,j,i)                &
2785                           + 5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
2786                   ENDDO
2787                ENDDO
2788             ENDDO
2789          ENDIF
2790       ENDIF
2791
2792    ENDDO
2793
2794    CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )
2795
2796 END SUBROUTINE chem_prognostic_equations
2797
2798
2799!------------------------------------------------------------------------------!
2800! Description:
2801! ------------
2802!> Subroutine calculating prognostic equations for chemical species
2803!> (cache-optimized).
2804!> Routine is called separately for each chemical species over a loop from
2805!> prognostic_equations.
2806!------------------------------------------------------------------------------!
2807 SUBROUTINE chem_prognostic_equations_ij( i, j, i_omp_start, tn )
2808
2809
2810    INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn
2811    INTEGER(iwp) :: ilsp
2812!
2813!-- local variables
2814
2815    INTEGER :: k
2816
2817    DO  ilsp = 1, nvar
2818!
2819!--    Tendency-terms for chem spcs.
2820       tend(:,j,i) = 0.0_wp
2821!
2822!--    Advection terms
2823       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2824          IF ( ws_scheme_sca )  THEN
2825             CALL advec_s_ws( cs_advc_flags_s,                                                     &
2826                              i,                                                                   &
2827                              j,                                                                   &
2828                              chem_species(ilsp)%conc,                                             &
2829                              'kc',                                                                &
2830                              chem_species(ilsp)%flux_s_cs,                                        &
2831                              chem_species(ilsp)%diss_s_cs,                                        &
2832                              chem_species(ilsp)%flux_l_cs,                                        &
2833                              chem_species(ilsp)%diss_l_cs,                                        &
2834                              i_omp_start,                                                         &
2835                              tn,                                                                  &
2836                              bc_dirichlet_l  .OR.  bc_radiation_l  .OR.  decycle_chem_lr,         &
2837                              bc_dirichlet_n  .OR.  bc_radiation_n  .OR.  decycle_chem_ns,         &
2838                              bc_dirichlet_r  .OR.  bc_radiation_r  .OR.  decycle_chem_lr,         &
2839                              bc_dirichlet_s  .OR.  bc_radiation_s  .OR.  decycle_chem_ns,         &
2840                              monotonic_limiter_z )
2841          ELSE
2842             CALL advec_s_pw( i, j, chem_species(ilsp)%conc )
2843          ENDIF
2844       ELSE
2845          CALL advec_s_up( i, j, chem_species(ilsp)%conc )
2846       ENDIF
2847!
2848!--    Diffusion terms (the last three arguments are zero)
2849
2850       CALL diffusion_s( i, j, chem_species(ilsp)%conc,                                            &
2851            surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),                              &
2852            surf_def_h(2)%cssws(ilsp,:),                                                           &
2853            surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:),                                    &
2854            surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),                              &
2855            surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),                              &
2856            surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:),                              &
2857            surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:),                              &
2858            surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),                              &
2859            surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
2860!
2861!--    Prognostic equation for chem spcs
2862       DO  k = nzb+1, nzt
2863          chem_species(ilsp)%conc_p(k,j,i) = chem_species(ilsp)%conc(k,j,i) + ( dt_3d  *           &
2864               ( tsc(2) * tend(k,j,i) +                                                            &
2865               tsc(3) * chem_species(ilsp)%tconc_m(k,j,i) )                                        &
2866               - tsc(5) * rdf_sc(k)                                                                &
2867               * ( chem_species(ilsp)%conc(k,j,i) - chem_species(ilsp)%conc_pr_init(k) )           &
2868               )                                                                                   &
2869               * MERGE( 1.0_wp, 0.0_wp,                                                            &
2870               BTEST( wall_flags_total_0(k,j,i), 0 )                                               &
2871               )
2872
2873          IF ( chem_species(ilsp)%conc_p(k,j,i) < 0.0_wp )  THEN               
2874             chem_species(ilsp)%conc_p(k,j,i) = 0.1_wp * chem_species(ilsp)%conc(k,j,i)    !FKS6               
2875          ENDIF
2876       ENDDO
2877!
2878!--    Calculate tendencies for the next Runge-Kutta step
2879       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2880          IF ( intermediate_timestep_count == 1 )  THEN
2881             DO  k = nzb+1, nzt
2882                chem_species(ilsp)%tconc_m(k,j,i) = tend(k,j,i)
2883             ENDDO
2884          ELSEIF ( intermediate_timestep_count < &
2885               intermediate_timestep_count_max )  THEN
2886             DO  k = nzb+1, nzt
2887                chem_species(ilsp)%tconc_m(k,j,i) = -9.5625_wp * tend(k,j,i) +                     &
2888                     5.3125_wp * chem_species(ilsp)%tconc_m(k,j,i)
2889             ENDDO
2890          ENDIF
2891       ENDIF
2892
2893    ENDDO
2894
2895 END SUBROUTINE chem_prognostic_equations_ij
2896
2897
2898!------------------------------------------------------------------------------!
2899! Description:
2900! ------------
2901!> Subroutine to read restart data of chemical species
2902!------------------------------------------------------------------------------!
2903 SUBROUTINE chem_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
2904                            nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,   &
2905                            nys_on_file, tmp_3d, found )
2906
2907    USE control_parameters
2908
2909
2910    CHARACTER (LEN=20) :: spc_name_av !<   
2911
2912    INTEGER(iwp) ::  lsp             !<
2913    INTEGER(iwp) ::  k               !<
2914    INTEGER(iwp) ::  nxlc            !<
2915    INTEGER(iwp) ::  nxlf            !<
2916    INTEGER(iwp) ::  nxl_on_file     !<   
2917    INTEGER(iwp) ::  nxrc            !<
2918    INTEGER(iwp) ::  nxrf            !<
2919    INTEGER(iwp) ::  nxr_on_file     !<   
2920    INTEGER(iwp) ::  nync            !<
2921    INTEGER(iwp) ::  nynf            !<
2922    INTEGER(iwp) ::  nyn_on_file     !<   
2923    INTEGER(iwp) ::  nysc            !<
2924    INTEGER(iwp) ::  nysf            !<
2925    INTEGER(iwp) ::  nys_on_file     !<   
2926
2927    LOGICAL, INTENT(OUT) :: found
2928
2929    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
2930
2931
2932    found = .FALSE. 
2933
2934
2935    IF ( ALLOCATED(chem_species) )  THEN
2936
2937       DO  lsp = 1, nspec
2938
2939          !< for time-averaged chemical conc.
2940          spc_name_av  =  TRIM( chem_species(lsp)%name )//'_av'
2941
2942          IF ( restart_string(1:length) == TRIM( chem_species(lsp)%name) )    &
2943             THEN
2944             !< read data into tmp_3d
2945             IF ( k == 1 )  READ ( 13 )  tmp_3d 
2946             !< fill ..%conc in the restart run   
2947             chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                    &
2948                  nxlc-nbgp:nxrc+nbgp) =                  & 
2949                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2950             found = .TRUE.
2951          ELSEIF (restart_string(1:length) == spc_name_av )  THEN
2952             IF ( k == 1 )  READ ( 13 )  tmp_3d
2953             chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,                 &
2954                  nxlc-nbgp:nxrc+nbgp) =               &
2955                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2956             found = .TRUE.
2957          ENDIF
2958
2959       ENDDO
2960
2961    ENDIF
2962
2963
2964 END SUBROUTINE chem_rrd_local
2965
2966
2967!-------------------------------------------------------------------------------!
2968!> Description:
2969!> Calculation of horizontally averaged profiles
2970!> This routine is called for every statistic region (sr) defined by the user,
2971!> but at least for the region "total domain" (sr=0).
2972!> quantities.
2973!-------------------------------------------------------------------------------!
2974 SUBROUTINE chem_statistics( mode, sr, tn )
2975
2976
2977    USE arrays_3d
2978
2979    USE statistics
2980
2981
2982    CHARACTER (LEN=*) ::  mode   !<
2983
2984    INTEGER(iwp) ::  i    !< running index on x-axis
2985    INTEGER(iwp) ::  j    !< running index on y-axis
2986    INTEGER(iwp) ::  k    !< vertical index counter
2987    INTEGER(iwp) ::  sr   !< statistical region
2988    INTEGER(iwp) ::  tn   !< thread number
2989    INTEGER(iwp) ::  lpr  !< running index chem spcs
2990
2991    IF ( mode == 'profiles' )  THEN
2992       !
2993!
2994!--    Sample on how to calculate horizontally averaged profiles of user-
2995!--    defined quantities. Each quantity is identified by the index
2996!--    "pr_palm+#" where "#" is an integer starting from 1. These
2997!--    user-profile-numbers must also be assigned to the respective strings
2998!--    given by data_output_pr_cs in routine user_check_data_output_pr.
2999!--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
3000!--                     w*pt*), dim-4 = statistical region.
3001
3002!$OMP DO
3003       DO  i = nxl, nxr
3004          DO  j = nys, nyn
3005             DO  k = nzb, nzt+1
3006                DO lpr = 1, cs_pr_count
3007
3008                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
3009                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
3010                        rmask(j,i,sr)  *                                   &
3011                        MERGE( 1.0_wp, 0.0_wp,                             &
3012                        BTEST( wall_flags_total_0(k,j,i), 22 ) )
3013                ENDDO
3014             ENDDO
3015          ENDDO
3016       ENDDO
3017    ELSEIF ( mode == 'time_series' )  THEN
3018!      @todo
3019    ENDIF
3020
3021 END SUBROUTINE chem_statistics
3022
3023
3024!------------------------------------------------------------------------------!
3025! Description:
3026! ------------
3027!> Subroutine for swapping of timelevels for chemical species
3028!> called out from subroutine swap_timelevel
3029!------------------------------------------------------------------------------!
3030
3031
3032 SUBROUTINE chem_swap_timelevel( level )
3033
3034
3035    INTEGER(iwp), INTENT(IN) ::  level
3036!
3037!-- local variables
3038    INTEGER(iwp)             ::  lsp
3039
3040
3041    IF ( level == 0 )  THEN
3042       DO  lsp=1, nvar                                       
3043          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
3044          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
3045       ENDDO
3046    ELSE
3047       DO  lsp=1, nvar                                       
3048          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
3049          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
3050       ENDDO
3051    ENDIF
3052
3053    RETURN
3054 END SUBROUTINE chem_swap_timelevel
3055
3056
3057!------------------------------------------------------------------------------!
3058! Description:
3059! ------------
3060!> Subroutine to write restart data for chemistry model
3061!------------------------------------------------------------------------------!
3062 SUBROUTINE chem_wrd_local
3063
3064
3065    INTEGER(iwp) ::  lsp  !< running index for chem spcs.
3066
3067    DO  lsp = 1, nspec
3068       CALL wrd_write_string( TRIM( chem_species(lsp)%name ) )
3069       WRITE ( 14 )  chem_species(lsp)%conc
3070       CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
3071       WRITE ( 14 )  chem_species(lsp)%conc_av
3072    ENDDO
3073
3074 END SUBROUTINE chem_wrd_local
3075
3076
3077!-------------------------------------------------------------------------------!
3078! Description:
3079! ------------
3080!> Subroutine to calculate the deposition of gases and PMs. For now deposition
3081!> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
3082!> considered. The deposition of particles is derived following Zhang et al.,
3083!> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
3084!>     
3085!> @TODO: Consider deposition on vertical surfaces   
3086!> @TODO: Consider overlaying horizontal surfaces
3087!> @TODO: Consider resolved vegetation   
3088!> @TODO: Check error messages
3089!-------------------------------------------------------------------------------!
3090 SUBROUTINE chem_depo( i, j )
3091
3092    USE control_parameters,                                                 &   
3093         ONLY:  dt_3d, intermediate_timestep_count, latitude,               &
3094                time_since_reference_point
3095
3096    USE arrays_3d,                                                          &
3097         ONLY:  dzw, rho_air_zw
3098
3099    USE palm_date_time_mod,                                                 &
3100         ONLY:  get_date_time
3101
3102    USE surface_mod,                                                        &
3103         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h,        &
3104         surf_type, surf_usm_h
3105
3106    USE radiation_model_mod,                                                &
3107         ONLY:  cos_zenith
3108
3109
3110    INTEGER(iwp) ::  day_of_year                   !< current day of the year
3111    INTEGER(iwp), INTENT(IN) ::  i
3112    INTEGER(iwp), INTENT(IN) ::  j
3113    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
3114    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
3115    INTEGER(iwp) ::  luv_palm                      !< index of PALM LSM vegetation_type at current surface element
3116    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
3117    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current surface element
3118    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current surface element
3119    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current surface element
3120    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface element
3121    INTEGER(iwp) ::  luv_dep                       !< matching DEPAC LU to luv_palm
3122    INTEGER(iwp) ::  lup_dep                       !< matching DEPAC LU to lup_palm
3123    INTEGER(iwp) ::  luw_dep                       !< matching DEPAC LU to luw_palm
3124    INTEGER(iwp) ::  luu_dep                       !< matching DEPAC LU to luu_palm
3125    INTEGER(iwp) ::  lug_dep                       !< matching DEPAC LU to lug_palm
3126    INTEGER(iwp) ::  lud_dep                       !< matching DEPAC LU to lud_palm
3127    INTEGER(iwp) ::  m                             !< index for horizontal surfaces
3128
3129    INTEGER(iwp) ::  pspec                         !< running index
3130    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
3131!
3132!-- Vegetation                                               !< Assign PALM classes to DEPAC land use classes
3133    INTEGER(iwp) ::  ind_luv_user = 0                        !<  ERROR as no class given in PALM
3134    INTEGER(iwp) ::  ind_luv_b_soil = 1                      !<  assigned to ilu_desert
3135    INTEGER(iwp) ::  ind_luv_mixed_crops = 2                 !<  assigned to ilu_arable
3136    INTEGER(iwp) ::  ind_luv_s_grass = 3                     !<  assigned to ilu_grass
3137    INTEGER(iwp) ::  ind_luv_ev_needle_trees = 4             !<  assigned to ilu_coniferous_forest
3138    INTEGER(iwp) ::  ind_luv_de_needle_trees = 5             !<  assigned to ilu_coniferous_forest
3139    INTEGER(iwp) ::  ind_luv_ev_broad_trees = 6              !<  assigned to ilu_tropical_forest
3140    INTEGER(iwp) ::  ind_luv_de_broad_trees = 7              !<  assigned to ilu_deciduous_forest
3141    INTEGER(iwp) ::  ind_luv_t_grass = 8                     !<  assigned to ilu_grass
3142    INTEGER(iwp) ::  ind_luv_desert = 9                      !<  assigned to ilu_desert
3143    INTEGER(iwp) ::  ind_luv_tundra = 10                     !<  assigned to ilu_other
3144    INTEGER(iwp) ::  ind_luv_irr_crops = 11                  !<  assigned to ilu_arable
3145    INTEGER(iwp) ::  ind_luv_semidesert = 12                 !<  assigned to ilu_other
3146    INTEGER(iwp) ::  ind_luv_ice = 13                        !<  assigned to ilu_ice
3147    INTEGER(iwp) ::  ind_luv_marsh = 14                      !<  assigned to ilu_other
3148    INTEGER(iwp) ::  ind_luv_ev_shrubs = 15                  !<  assigned to ilu_mediterrean_scrub
3149    INTEGER(iwp) ::  ind_luv_de_shrubs = 16                  !<  assigned to ilu_mediterrean_scrub
3150    INTEGER(iwp) ::  ind_luv_mixed_forest = 17               !<  assigned to ilu_coniferous_forest (ave(decid+conif))
3151    INTEGER(iwp) ::  ind_luv_intrup_forest = 18              !<  assigned to ilu_other (ave(other+decid))
3152!
3153!-- Water
3154    INTEGER(iwp) ::  ind_luw_user = 0                        !<  ERROR as no class given in PALM 
3155    INTEGER(iwp) ::  ind_luw_lake = 1                        !<  assigned to ilu_water_inland
3156    INTEGER(iwp) ::  ind_luw_river = 2                       !<  assigned to ilu_water_inland
3157    INTEGER(iwp) ::  ind_luw_ocean = 3                       !<  assigned to ilu_water_sea
3158    INTEGER(iwp) ::  ind_luw_pond = 4                        !<  assigned to ilu_water_inland
3159    INTEGER(iwp) ::  ind_luw_fountain = 5                    !<  assigned to ilu_water_inland
3160!
3161!-- Pavement
3162    INTEGER(iwp) ::  ind_lup_user = 0                        !<  ERROR as no class given in PALM
3163    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !<  assigned to ilu_desert
3164    INTEGER(iwp) ::  ind_lup_asph = 2                        !<  assigned to ilu_desert
3165    INTEGER(iwp) ::  ind_lup_conc = 3                        !<  assigned to ilu_desert
3166    INTEGER(iwp) ::  ind_lup_sett = 4                        !<  assigned to ilu_desert
3167    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !<  assigned to ilu_desert
3168    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !<  assigned to ilu_desert
3169    INTEGER(iwp) ::  ind_lup_metal = 7                       !<  assigned to ilu_desert
3170    INTEGER(iwp) ::  ind_lup_wood = 8                        !<  assigned to ilu_desert
3171    INTEGER(iwp) ::  ind_lup_gravel = 9                      !<  assigned to ilu_desert
3172    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !<  assigned to ilu_desert
3173    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !<  assigned to ilu_desert
3174    INTEGER(iwp) ::  ind_lup_woodchips = 12                  !<  assigned to ilu_desert
3175    INTEGER(iwp) ::  ind_lup_tartan = 13                     !<  assigned to ilu_desert
3176    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !<  assigned to ilu_desert
3177    INTEGER(iwp) ::  ind_lup_clay = 15                       !<  assigned to ilu_desert
3178!
3179!-- Particle parameters according to the respective aerosol classes (PM25, PM10)
3180    INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
3181    INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
3182    INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
3183
3184    INTEGER(iwp) ::  part_type          !< index for particle type (PM10 or PM25) in particle_pars
3185
3186    INTEGER(iwp) ::  nwet               !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
3187
3188    REAL(wp) ::  dt_chem                !< length of chem time step
3189    REAL(wp) ::  dh                     !< vertical grid size
3190    REAL(wp) ::  inv_dh                 !< inverse of vertical grid size
3191    REAL(wp) ::  dt_dh                  !< dt_chem/dh
3192
3193    REAL(wp) ::  dens              !< density at layer k at i,j 
3194    REAL(wp) ::  r_aero_surf       !< aerodynamic resistance (s/m) at current surface element
3195    REAL(wp) ::  ustar_surf        !< ustar at current surface element
3196    REAL(wp) ::  z0h_surf          !< roughness length for heat at current surface element
3197    REAL(wp) ::  solar_rad         !< solar radiation, direct and diffuse, at current surface element
3198    REAL(wp) ::  ppm2ugm3          !< conversion factor from ppm to ug/m3
3199    REAL(wp) ::  rh_surf           !< relative humidity at current surface element
3200    REAL(wp) ::  lai               !< leaf area index at current surface element
3201    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
3202
3203    REAL(wp) ::  slinnfac       
3204    REAL(wp) ::  visc              !< Viscosity
3205    REAL(wp) ::  vs                !< Sedimentation velocity
3206    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
3207    REAL(wp) ::  rs                !< Sedimentaion resistance (s/m)
3208    REAL(wp) ::  rb                !< quasi-laminar boundary layer resistance (s/m)
3209    REAL(wp) ::  rc_tot            !< total canopy resistance (s/m)
3210
3211    REAL(wp) ::  conc_ijk_ugm3     !< concentration at i, j, k in ug/m3
3212    REAL(wp) ::  diffusivity       !< diffusivity
3213
3214
3215    REAL(wp), DIMENSION(nspec) ::  bud_luv      !< budget for LSM vegetation type at current surface element
3216    REAL(wp), DIMENSION(nspec) ::  bud_lup      !< budget for LSM pavement type at current surface element
3217    REAL(wp), DIMENSION(nspec) ::  bud_luw      !< budget for LSM water type at current surface element
3218    REAL(wp), DIMENSION(nspec) ::  bud_luu      !< budget for USM walls/roofs at current surface element
3219    REAL(wp), DIMENSION(nspec) ::  bud_lug      !< budget for USM green surfaces at current surface element
3220    REAL(wp), DIMENSION(nspec) ::  bud_lud      !< budget for USM windows at current surface element
3221    REAL(wp), DIMENSION(nspec) ::  bud          !< overall budget at current surface element
3222    REAL(wp), DIMENSION(nspec) ::  conc_ijk     !< concentration at i,j,k
3223    REAL(wp), DIMENSION(nspec) ::  ccomp_tot    !< total compensation point (ug/m3), for now kept to zero for all species!
3224                                                 
3225
3226    REAL(wp) ::  temp_tmp       !< temperatur at i,j,k
3227    REAL(wp) ::  ts             !< surface temperatur in degrees celsius
3228    REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
3229!
3230!-- Particle parameters (PM10 (1), PM25 (2))
3231!-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
3232!-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
3233    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/ &
3234         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &  !<  1
3235         0.7e-6_wp, 1.14e3_wp, 1.082_wp &   !<  2
3236         /), (/ 3, 2 /) )
3237
3238    LOGICAL ::  match_lsm     !< flag indicating natural-type surface
3239    LOGICAL ::  match_usm     !< flag indicating urban-type surface
3240!
3241!-- List of names of possible tracers
3242    CHARACTER(LEN=*), PARAMETER ::  pspecnames(nposp) = (/ &
3243         'NO2           ', &    !< NO2
3244         'NO            ', &    !< NO
3245         'O3            ', &    !< O3
3246         'CO            ', &    !< CO
3247         'form          ', &    !< FORM
3248         'ald           ', &    !< ALD
3249         'pan           ', &    !< PAN
3250         'mgly          ', &    !< MGLY
3251         'par           ', &    !< PAR
3252         'ole           ', &    !< OLE
3253         'eth           ', &    !< ETH
3254         'tol           ', &    !< TOL
3255         'cres          ', &    !< CRES
3256         'xyl           ', &    !< XYL
3257         'SO4a_f        ', &    !< SO4a_f
3258         'SO2           ', &    !< SO2
3259         'HNO2          ', &    !< HNO2
3260         'CH4           ', &    !< CH4
3261         'NH3           ', &    !< NH3
3262         'NO3           ', &    !< NO3
3263         'OH            ', &    !< OH
3264         'HO2           ', &    !< HO2
3265         'N2O5          ', &    !< N2O5
3266         'SO4a_c        ', &    !< SO4a_c
3267         'NH4a_f        ', &    !< NH4a_f
3268         'NO3a_f        ', &    !< NO3a_f
3269         'NO3a_c        ', &    !< NO3a_c
3270         'C2O3          ', &    !< C2O3
3271         'XO2           ', &    !< XO2
3272         'XO2N          ', &    !< XO2N
3273         'cro           ', &    !< CRO
3274         'HNO3          ', &    !< HNO3
3275         'H2O2          ', &    !< H2O2
3276         'iso           ', &    !< ISO
3277         'ispd          ', &    !< ISPD
3278         'to2           ', &    !< TO2
3279         'open          ', &    !< OPEN
3280         'terp          ', &    !< TERP
3281         'ec_f          ', &    !< EC_f
3282         'ec_c          ', &    !< EC_c
3283         'pom_f         ', &    !< POM_f
3284         'pom_c         ', &    !< POM_c
3285         'ppm_f         ', &    !< PPM_f
3286         'ppm_c         ', &    !< PPM_c
3287         'na_ff         ', &    !< Na_ff
3288         'na_f          ', &    !< Na_f
3289         'na_c          ', &    !< Na_c
3290         'na_cc         ', &    !< Na_cc
3291         'na_ccc        ', &    !< Na_ccc
3292         'dust_ff       ', &    !< dust_ff
3293         'dust_f        ', &    !< dust_f
3294         'dust_c        ', &    !< dust_c
3295         'dust_cc       ', &    !< dust_cc
3296         'dust_ccc      ', &    !< dust_ccc
3297         'tpm10         ', &    !< tpm10
3298         'tpm25         ', &    !< tpm25
3299         'tss           ', &    !< tss
3300         'tdust         ', &    !< tdust
3301         'tc            ', &    !< tc
3302         'tcg           ', &    !< tcg
3303         'tsoa          ', &    !< tsoa
3304         'tnmvoc        ', &    !< tnmvoc
3305         'SOxa          ', &    !< SOxa
3306         'NOya          ', &    !< NOya
3307         'NHxa          ', &    !< NHxa
3308         'NO2_obs       ', &    !< NO2_obs
3309         'tpm10_biascorr', &    !< tpm10_biascorr
3310         'tpm25_biascorr', &    !< tpm25_biascorr
3311         'O3_biascorr   ' /)    !< o3_biascorr
3312!
3313!-- tracer mole mass:
3314    REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
3315         xm_O * 2 + xm_N, &                         !< NO2
3316         xm_O + xm_N, &                             !< NO
3317         xm_O * 3, &                                !< O3
3318         xm_C + xm_O, &                             !< CO
3319         xm_H * 2 + xm_C + xm_O, &                  !< FORM
3320         xm_H * 3 + xm_C * 2 + xm_O, &              !< ALD
3321         xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, &   !< PAN
3322         xm_H * 4 + xm_C * 3 + xm_O * 2, &          !< MGLY
3323         xm_H * 3 + xm_C, &                         !< PAR
3324         xm_H * 3 + xm_C * 2, &                     !< OLE
3325         xm_H * 4 + xm_C * 2, &                     !< ETH
3326         xm_H * 8 + xm_C * 7, &                     !< TOL
3327         xm_H * 8 + xm_C * 7 + xm_O, &              !< CRES
3328         xm_H * 10 + xm_C * 8, &                    !< XYL
3329         xm_S + xm_O * 4, &                         !< SO4a_f
3330         xm_S + xm_O * 2, &                         !< SO2
3331         xm_H + xm_O * 2 + xm_N, &                  !< HNO2
3332         xm_H * 4 + xm_C, &                         !< CH4
3333         xm_H * 3 + xm_N, &                         !< NH3
3334         xm_O * 3 + xm_N, &                         !< NO3
3335         xm_H + xm_O, &                             !< OH
3336         xm_H + xm_O * 2, &                         !< HO2
3337         xm_O * 5 + xm_N * 2, &                     !< N2O5
3338         xm_S + xm_O * 4, &                         !< SO4a_c
3339         xm_H * 4 + xm_N, &                         !< NH4a_f
3340         xm_O * 3 + xm_N, &                         !< NO3a_f
3341         xm_O * 3 + xm_N, &                         !< NO3a_c
3342         xm_C * 2 + xm_O * 3, &                     !< C2O3
3343         xm_dummy, &                                !< XO2
3344         xm_dummy, &                                !< XO2N
3345         xm_dummy, &                                !< CRO
3346         xm_H + xm_O * 3 + xm_N, &                  !< HNO3
3347         xm_H * 2 + xm_O * 2, &                     !< H2O2
3348         xm_H * 8 + xm_C * 5, &                     !< ISO
3349         xm_dummy, &                                !< ISPD
3350         xm_dummy, &                                !< TO2
3351         xm_dummy, &                                !< OPEN
3352         xm_H * 16 + xm_C * 10, &                   !< TERP
3353         xm_dummy, &                                !< EC_f
3354         xm_dummy, &                                !< EC_c
3355         xm_dummy, &                                !< POM_f
3356         xm_dummy, &                                !< POM_c
3357         xm_dummy, &                                !< PPM_f
3358         xm_dummy, &                                !< PPM_c
3359         xm_Na, &                                   !< Na_ff
3360         xm_Na, &                                   !< Na_f
3361         xm_Na, &                                   !< Na_c
3362         xm_Na, &                                   !< Na_cc
3363         xm_Na, &                                   !< Na_ccc
3364         xm_dummy, &                                !< dust_ff
3365         xm_dummy, &                                !< dust_f
3366         xm_dummy, &                                !< dust_c
3367         xm_dummy, &                                !< dust_cc
3368         xm_dummy, &                                !< dust_ccc
3369         xm_dummy, &                                !< tpm10
3370         xm_dummy, &                                !< tpm25
3371         xm_dummy, &                                !< tss
3372         xm_dummy, &                                !< tdust
3373         xm_dummy, &                                !< tc
3374         xm_dummy, &                                !< tcg
3375         xm_dummy, &                                !< tsoa
3376         xm_dummy, &                                !< tnmvoc
3377         xm_dummy, &                                !< SOxa
3378         xm_dummy, &                                !< NOya
3379         xm_dummy, &                                !< NHxa
3380         xm_O * 2 + xm_N, &                         !< NO2_obs
3381         xm_dummy, &                                !< tpm10_biascorr
3382         xm_dummy, &                                !< tpm25_biascorr
3383         xm_O * 3 /)                                !< o3_biascorr
3384!
3385!-- Get current day of the year
3386    CALL get_date_time( time_since_reference_point, day_of_year = day_of_year )
3387!
3388!-- Initialize surface element m
3389    m = 0
3390    k = 0
3391!
3392!-- LSM or USM surface present at i,j:
3393!-- Default surfaces are NOT considered for deposition
3394    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
3395    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
3396!
3397!--For LSM surfaces
3398
3399    IF ( match_lsm )  THEN
3400!
3401!--    Get surface element information at i,j:
3402       m = surf_lsm_h%start_index(j,i)
3403       k = surf_lsm_h%k(m)
3404!
3405!--    Get needed variables for surface element m
3406       ustar_surf  = surf_lsm_h%us(m)
3407       z0h_surf    = surf_lsm_h%z0h(m)
3408       r_aero_surf = surf_lsm_h%r_a(m)
3409       solar_rad   = surf_lsm_h%rad_sw_dir(m) + surf_lsm_h%rad_sw_dif(m)
3410       
3411       lai = surf_lsm_h%lai(m)
3412       sai = lai + 1
3413!
3414!--    For small grid spacing neglect R_a
3415       IF ( dzw(k) <= 1.0 )  THEN
3416          r_aero_surf = 0.0_wp
3417       ENDIF
3418!
3419!--    Initialize lu's
3420       luv_palm = 0
3421       luv_dep = 0
3422       lup_palm = 0
3423       lup_dep = 0
3424       luw_palm = 0
3425       luw_dep = 0
3426!
3427!--    Initialize budgets
3428       bud_luv = 0.0_wp
3429       bud_lup = 0.0_wp
3430       bud_luw = 0.0_wp
3431!
3432!--    Get land use for i,j and assign to DEPAC lu
3433       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
3434          luv_palm = surf_lsm_h%vegetation_type(m)
3435          IF ( luv_palm == ind_luv_user )  THEN
3436             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
3437             CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
3438          ELSEIF ( luv_palm == ind_luv_b_soil )  THEN
3439             luv_dep = 9
3440          ELSEIF ( luv_palm == ind_luv_mixed_crops )  THEN
3441             luv_dep = 2
3442          ELSEIF ( luv_palm == ind_luv_s_grass )  THEN
3443             luv_dep = 1
3444          ELSEIF ( luv_palm == ind_luv_ev_needle_trees )  THEN
3445             luv_dep = 4
3446          ELSEIF ( luv_palm == ind_luv_de_needle_trees )  THEN
3447             luv_dep = 4
3448          ELSEIF ( luv_palm == ind_luv_ev_broad_trees )  THEN
3449             luv_dep = 12
3450          ELSEIF ( luv_palm == ind_luv_de_broad_trees )  THEN
3451             luv_dep = 5
3452          ELSEIF ( luv_palm == ind_luv_t_grass )  THEN
3453             luv_dep = 1
3454          ELSEIF ( luv_palm == ind_luv_desert )  THEN
3455             luv_dep = 9
3456          ELSEIF ( luv_palm == ind_luv_tundra )  THEN
3457             luv_dep = 8
3458          ELSEIF ( luv_palm == ind_luv_irr_crops )  THEN
3459             luv_dep = 2
3460          ELSEIF ( luv_palm == ind_luv_semidesert )  THEN
3461             luv_dep = 8
3462          ELSEIF ( luv_palm == ind_luv_ice )  THEN
3463             luv_dep = 10
3464          ELSEIF ( luv_palm == ind_luv_marsh )  THEN
3465             luv_dep = 8
3466          ELSEIF ( luv_palm == ind_luv_ev_shrubs )  THEN
3467             luv_dep = 14
3468          ELSEIF ( luv_palm == ind_luv_de_shrubs )  THEN
3469             luv_dep = 14
3470          ELSEIF ( luv_palm == ind_luv_mixed_forest )  THEN
3471             luv_dep = 4
3472          ELSEIF ( luv_palm == ind_luv_intrup_forest )  THEN
3473             luv_dep = 8     
3474          ENDIF
3475       ENDIF
3476
3477       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
3478          lup_palm = surf_lsm_h%pavement_type(m)
3479          IF ( lup_palm == ind_lup_user )  THEN
3480             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
3481             CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 )
3482          ELSEIF ( lup_palm == ind_lup_asph_conc )  THEN
3483             lup_dep = 9
3484          ELSEIF ( lup_palm == ind_lup_asph )  THEN
3485             lup_dep = 9
3486          ELSEIF ( lup_palm == ind_lup_conc )  THEN
3487             lup_dep = 9
3488          ELSEIF ( lup_palm == ind_lup_sett )  THEN
3489             lup_dep = 9
3490          ELSEIF ( lup_palm == ind_lup_pav_stones )  THEN
3491             lup_dep = 9
3492          ELSEIF ( lup_palm == ind_lup_cobblest )  THEN
3493             lup_dep = 9       
3494          ELSEIF ( lup_palm == ind_lup_metal )  THEN
3495             lup_dep = 9
3496          ELSEIF ( lup_palm == ind_lup_wood )  THEN
3497             lup_dep = 9   
3498          ELSEIF ( lup_palm == ind_lup_gravel )  THEN
3499             lup_dep = 9
3500          ELSEIF ( lup_palm == ind_lup_f_gravel )  THEN
3501             lup_dep = 9
3502          ELSEIF ( lup_palm == ind_lup_pebblest )  THEN
3503             lup_dep = 9
3504          ELSEIF ( lup_palm == ind_lup_woodchips )  THEN
3505             lup_dep = 9
3506          ELSEIF ( lup_palm == ind_lup_tartan )  THEN
3507             lup_dep = 9
3508          ELSEIF ( lup_palm == ind_lup_art_turf )  THEN
3509             lup_dep = 9
3510          ELSEIF ( lup_palm == ind_lup_clay )  THEN
3511             lup_dep = 9
3512          ENDIF
3513       ENDIF
3514
3515       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
3516          luw_palm = surf_lsm_h%water_type(m)     
3517          IF ( luw_palm == ind_luw_user )  THEN
3518             message_string = 'No water type defined. Please define water type to enable deposition calculation'
3519             CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 )
3520          ELSEIF ( luw_palm ==  ind_luw_lake )  THEN
3521             luw_dep = 13
3522          ELSEIF ( luw_palm == ind_luw_river )  THEN
3523             luw_dep = 13
3524          ELSEIF ( luw_palm == ind_luw_ocean )  THEN
3525             luw_dep = 6
3526          ELSEIF ( luw_palm == ind_luw_pond )  THEN
3527             luw_dep = 13
3528          ELSEIF ( luw_palm == ind_luw_fountain )  THEN
3529             luw_dep = 13 
3530          ENDIF
3531       ENDIF
3532!
3533!--    Set wetness indicator to dry or wet for lsm vegetation or pavement
3534       IF ( surf_lsm_h%c_liq(m) > 0 )  THEN
3535          nwet = 1
3536       ELSE
3537          nwet = 0
3538       ENDIF
3539!
3540!--    Compute length of time step
3541       IF ( call_chem_at_all_substeps )  THEN
3542          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
3543       ELSE
3544          dt_chem = dt_3d
3545       ENDIF
3546
3547       dh = dzw(k)
3548       inv_dh = 1.0_wp / dh
3549       dt_dh = dt_chem/dh
3550!
3551!--    Concentration at i,j,k
3552       DO  lsp = 1, nspec
3553          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
3554       ENDDO
3555
3556!--    Temperature at i,j,k
3557       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
3558
3559       ts       = temp_tmp - 273.15  !< in degrees celcius
3560!
3561!--    Viscosity of air
3562       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
3563!
3564!--    Air density at k
3565       dens = rho_air_zw(k)
3566!
3567!--    Calculate relative humidity from specific humidity for DEPAC
3568       qv_tmp = MAX(q(k,j,i),0.0_wp)
3569       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
3570!
3571!-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
3572!-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
3573!
3574!--    Vegetation
3575       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
3576
3577!
3578!--       No vegetation on bare soil, desert or ice:
3579          IF ( ( luv_palm == ind_luv_b_soil ) .OR. &
3580                ( luv_palm == ind_luv_desert ) .OR. &
3581                 ( luv_palm == ind_luv_ice ) ) THEN
3582
3583             lai = 0.0_wp
3584             sai = 0.0_wp
3585
3586          ENDIF
3587         
3588          slinnfac = 1.0_wp
3589!
3590!--       Get deposition velocity vd
3591          DO  lsp = 1, nvar
3592!
3593!--          Initialize
3594             vs = 0.0_wp
3595             vd_lu = 0.0_wp
3596             rs = 0.0_wp
3597             rb = 0.0_wp
3598             rc_tot = 0.0_wp
3599             IF ( spc_names(lsp) == 'PM10' )  THEN
3600                part_type = 1
3601!
3602!--             Sedimentation velocity
3603                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3604                     particle_pars(ind_p_size, part_type), &
3605                     particle_pars(ind_p_slip, part_type), &
3606                     visc)
3607
3608                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3609                     vs, &
3610                     particle_pars(ind_p_size, part_type), &
3611                     particle_pars(ind_p_slip, part_type), &
3612                     nwet, temp_tmp, dens, visc, &
3613                     luv_dep,  &
3614                     r_aero_surf, ustar_surf )
3615
3616                bud_luv(lsp) = - conc_ijk(lsp) * &
3617                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3618
3619
3620             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3621                part_type = 2
3622!
3623!--             Sedimentation velocity
3624                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3625                     particle_pars(ind_p_size, part_type), &
3626                     particle_pars(ind_p_slip, part_type), &
3627                     visc)
3628
3629                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3630                     vs, &
3631                     particle_pars(ind_p_size, part_type), &
3632                     particle_pars(ind_p_slip, part_type), &
3633                     nwet, temp_tmp, dens, visc, &
3634                     luv_dep , &
3635                     r_aero_surf, ustar_surf )
3636
3637                bud_luv(lsp) = - conc_ijk(lsp) * &
3638                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3639
3640             ELSE !< GASES
3641!
3642!--             Read spc_name of current species for gas parameter
3643                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
3644                   i_pspec = 0
3645                   DO  pspec = 1, nposp
3646                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3647                         i_pspec = pspec
3648                      END IF
3649                   ENDDO
3650
3651                ELSE
3652!
3653!--             For now species not deposited
3654                   CYCLE
3655                ENDIF
3656!
3657!--             Factor used for conversion from ppb to ug/m3 :
3658!--             ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3659!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3660!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
3661!--             thus:
3662!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3663!--             Use density at k:
3664
3665                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3666!
3667!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3668                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3669                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
3670!
3671!--             Diffusivity for DEPAC relevant gases
3672!--             Use default value
3673                diffusivity            = 0.11e-4
3674!
3675!--             overwrite with known coefficients of diffusivity from Massman (1998)
3676                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3677                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3678                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3679                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3680                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3681                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3682                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3683!
3684!--             Get quasi-laminar boundary layer resistance rb:
3685                CALL get_rb_cell( (luv_dep == ilu_water_sea) .OR. (luv_dep == ilu_water_inland), &
3686                     z0h_surf, ustar_surf, diffusivity, &
3687                     rb )
3688!
3689!--             Get rc_tot
3690                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf, solar_rad, cos_zenith, &
3691                     rh_surf, lai, sai, nwet, luv_dep, 2, rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity, &
3692                     r_aero_surf , rb )
3693!
3694!--             Calculate budget
3695                IF ( rc_tot <= 0.0 )  THEN
3696
3697                   bud_luv(lsp) = 0.0_wp
3698
3699                ELSE
3700
3701                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3702
3703                   bud_luv(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3704                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3705                ENDIF
3706
3707             ENDIF
3708          ENDDO
3709       ENDIF
3710!
3711!--    Pavement
3712       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
3713!
3714!--       No vegetation on pavements:
3715          lai = 0.0_wp
3716          sai = 0.0_wp
3717
3718          slinnfac = 1.0_wp
3719!
3720!--       Get vd
3721          DO  lsp = 1, nvar
3722!
3723!--       Initialize
3724             vs = 0.0_wp
3725             vd_lu = 0.0_wp
3726             rs = 0.0_wp
3727             rb = 0.0_wp
3728             rc_tot = 0.0_wp
3729             IF ( spc_names(lsp) == 'PM10' )  THEN
3730                part_type = 1
3731!
3732!--             Sedimentation velocity:
3733                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3734                     particle_pars(ind_p_size, part_type), &
3735                     particle_pars(ind_p_slip, part_type), &
3736                     visc)
3737
3738                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3739                     vs, &
3740                     particle_pars(ind_p_size, part_type), &
3741                     particle_pars(ind_p_slip, part_type), &
3742                     nwet, temp_tmp, dens, visc, &
3743                     lup_dep,  &
3744                     r_aero_surf, ustar_surf )
3745
3746                bud_lup(lsp) = - conc_ijk(lsp) * &
3747                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3748
3749
3750             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3751                part_type = 2
3752!
3753!--             Sedimentation velocity:
3754                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3755                     particle_pars(ind_p_size, part_type), &
3756                     particle_pars(ind_p_slip, part_type), &
3757                     visc)
3758
3759                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3760                     vs, &
3761                     particle_pars(ind_p_size, part_type), &
3762                     particle_pars(ind_p_slip, part_type), &
3763                     nwet, temp_tmp, dens, visc, &
3764                     lup_dep, &
3765                     r_aero_surf, ustar_surf )
3766
3767                bud_lup(lsp) = - conc_ijk(lsp) * &
3768                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3769
3770             ELSE  !<GASES
3771!
3772!--             Read spc_name of current species for gas parameter
3773
3774                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
3775                   i_pspec = 0
3776                   DO  pspec = 1, nposp
3777                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3778                         i_pspec = pspec
3779                      END IF
3780                   ENDDO
3781
3782                ELSE
3783!
3784!--                For now species not deposited
3785                   CYCLE
3786                ENDIF
3787!
3788!--             Factor used for conversion from ppb to ug/m3 :
3789!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3790!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3791!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
3792!--             thus:
3793!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3794!--             Use density at lowest layer:
3795
3796                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3797!
3798!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3799                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3800                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
3801!
3802!--             Diffusivity for DEPAC relevant gases
3803!--             Use default value
3804                diffusivity            = 0.11e-4
3805!
3806!--             overwrite with known coefficients of diffusivity from Massman (1998)
3807                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3808                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3809                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3810                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3811                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3812                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3813                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3814!
3815!--             Get quasi-laminar boundary layer resistance rb:
3816                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland),   &
3817                     z0h_surf, ustar_surf, diffusivity, rb )
3818!
3819!--             Get rc_tot
3820                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,      &
3821                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lup_dep, 2,    &
3822                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,              &
3823                                         r_aero_surf , rb )
3824!
3825!--             Calculate budget
3826                IF ( rc_tot <= 0.0 )  THEN
3827                   bud_lup(lsp) = 0.0_wp
3828                ELSE
3829                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3830                   bud_lup(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3831                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3832                ENDIF
3833
3834             ENDIF
3835          ENDDO
3836       ENDIF
3837!
3838!--    Water
3839       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
3840!
3841!--       No vegetation on water:
3842          lai = 0.0_wp
3843          sai = 0.0_wp
3844          slinnfac = 1.0_wp
3845!
3846!--       Get vd
3847          DO  lsp = 1, nvar
3848!
3849!--          Initialize
3850             vs = 0.0_wp
3851             vd_lu = 0.0_wp
3852             rs = 0.0_wp
3853             rb = 0.0_wp
3854             rc_tot = 0.0_wp 
3855             IF ( spc_names(lsp) == 'PM10' )  THEN
3856                part_type = 1
3857!
3858!--             Sedimentation velocity:
3859                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3860                     particle_pars(ind_p_size, part_type), &
3861                     particle_pars(ind_p_slip, part_type), &
3862                     visc)
3863
3864                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3865                     vs, &
3866                     particle_pars(ind_p_size, part_type), &
3867                     particle_pars(ind_p_slip, part_type), &
3868                     nwet, temp_tmp, dens, visc, &
3869                     luw_dep, &
3870                     r_aero_surf, ustar_surf )
3871
3872                bud_luw(lsp) = - conc_ijk(lsp) * &
3873                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3874
3875             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
3876                part_type = 2
3877!
3878!--             Sedimentation velocity:
3879                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
3880                     particle_pars(ind_p_size, part_type), &
3881                     particle_pars(ind_p_slip, part_type), &
3882                     visc)
3883
3884                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
3885                     vs, &
3886                     particle_pars(ind_p_size, part_type), &
3887                     particle_pars(ind_p_slip, part_type), &
3888                     nwet, temp_tmp, dens, visc, &
3889                     luw_dep, &
3890                     r_aero_surf, ustar_surf )
3891
3892                bud_luw(lsp) = - conc_ijk(lsp) * &
3893                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3894
3895             ELSE  !<GASES
3896!
3897!--             Read spc_name of current species for gas PARAMETER
3898
3899                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
3900                   i_pspec = 0
3901                   DO  pspec = 1, nposp
3902                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
3903                         i_pspec = pspec
3904                      END IF
3905                   ENDDO
3906
3907                ELSE
3908!
3909!--                For now species not deposited
3910                   CYCLE
3911                ENDIF
3912!
3913!--             Factor used for conversion from ppb to ug/m3 :
3914!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
3915!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
3916!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
3917!--             thus:
3918!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
3919!--             Use density at lowest layer:
3920
3921                ppm2ugm3 = (dens/xm_air) * 0.001_wp  !< (mole air)/m3
3922!
3923!--             Atmospheric concentration in DEPAC is requested in ug/m3:
3924!--                 ug/m3        ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
3925                conc_ijk_ugm3 = conc_ijk(lsp) * ppm2ugm3 *  specmolm(i_pspec)  ! in ug/m3
3926!
3927!--             Diffusivity for DEPAC relevant gases
3928!--             Use default value
3929                diffusivity            = 0.11e-4
3930!
3931!--             overwrite with known coefficients of diffusivity from Massman (1998)
3932                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
3933                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
3934                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
3935                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
3936                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
3937                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
3938                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
3939!
3940!--             Get quasi-laminar boundary layer resistance rb:
3941                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland),  &
3942                     z0h_surf, ustar_surf, diffusivity, rb )
3943
3944!--             Get rc_tot
3945                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         & 
3946                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luw_dep, 2,    &
3947                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
3948                                         r_aero_surf , rb )
3949!
3950!--             Calculate budget
3951                IF ( rc_tot <= 0.0 )  THEN
3952
3953                   bud_luw(lsp) = 0.0_wp
3954
3955                ELSE
3956
3957                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
3958
3959                   bud_luw(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
3960                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
3961                ENDIF
3962
3963             ENDIF
3964          ENDDO
3965       ENDIF
3966
3967
3968       bud = 0.0_wp
3969!
3970!--    Calculate overall budget for surface m and adapt concentration
3971       DO  lsp = 1, nspec
3972
3973          bud(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_luv(lsp) + &
3974               surf_lsm_h%frac(ind_pav_green,m) * bud_lup(lsp) + &
3975               surf_lsm_h%frac(ind_wat_win,m) * bud_luw(lsp)
3976!
3977!--       Compute new concentration:
3978          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
3979
3980          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, conc_ijk(lsp))
3981
3982       ENDDO
3983
3984    ENDIF
3985!
3986!-- For USM surfaces   
3987
3988    IF ( match_usm )  THEN
3989!
3990!--    Get surface element information at i,j:
3991       m = surf_usm_h%start_index(j,i)
3992       k = surf_usm_h%k(m)
3993!
3994!--    Get needed variables for surface element m
3995       ustar_surf  = surf_usm_h%us(m)
3996       z0h_surf    = surf_usm_h%z0h(m)
3997       r_aero_surf = surf_usm_h%r_a(m)
3998       solar_rad   = surf_usm_h%rad_sw_dir(m) + surf_usm_h%rad_sw_dif(m)
3999       lai = surf_usm_h%lai(m)
4000       sai = lai + 1
4001!
4002!--    For small grid spacing neglect R_a
4003       IF ( dzw(k) <= 1.0 )  THEN
4004          r_aero_surf = 0.0_wp
4005       ENDIF
4006!
4007!--    Initialize lu's
4008       luu_palm = 0
4009       luu_dep = 0
4010       lug_palm = 0
4011       lug_dep = 0
4012       lud_palm = 0
4013       lud_dep = 0
4014!
4015!--    Initialize budgets
4016       bud_luu  = 0.0_wp
4017       bud_lug = 0.0_wp
4018       bud_lud = 0.0_wp
4019!
4020!--    Get land use for i,j and assign to DEPAC lu
4021       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
4022!
4023!--       For green urban surfaces (e.g. green roofs
4024!--       assume LU short grass
4025          lug_palm = ind_luv_s_grass
4026          IF ( lug_palm == ind_luv_user )  THEN
4027             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
4028             CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
4029          ELSEIF ( lug_palm == ind_luv_b_soil )  THEN
4030             lug_dep = 9
4031          ELSEIF ( lug_palm == ind_luv_mixed_crops )  THEN
4032             lug_dep = 2
4033          ELSEIF ( lug_palm == ind_luv_s_grass )  THEN
4034             lug_dep = 1
4035          ELSEIF ( lug_palm == ind_luv_ev_needle_trees )  THEN
4036             lug_dep = 4
4037          ELSEIF ( lug_palm == ind_luv_de_needle_trees )  THEN
4038             lug_dep = 4
4039          ELSEIF ( lug_palm == ind_luv_ev_broad_trees )  THEN
4040             lug_dep = 12
4041          ELSEIF ( lug_palm == ind_luv_de_broad_trees )  THEN
4042             lug_dep = 5
4043          ELSEIF ( lug_palm == ind_luv_t_grass )  THEN
4044             lug_dep = 1
4045          ELSEIF ( lug_palm == ind_luv_desert )  THEN
4046             lug_dep = 9
4047          ELSEIF ( lug_palm == ind_luv_tundra  )  THEN
4048             lug_dep = 8
4049          ELSEIF ( lug_palm == ind_luv_irr_crops )  THEN
4050             lug_dep = 2
4051          ELSEIF ( lug_palm == ind_luv_semidesert )  THEN
4052             lug_dep = 8
4053          ELSEIF ( lug_palm == ind_luv_ice )  THEN
4054             lug_dep = 10
4055          ELSEIF ( lug_palm == ind_luv_marsh )  THEN
4056             lug_dep = 8
4057          ELSEIF ( lug_palm == ind_luv_ev_shrubs )  THEN
4058             lug_dep = 14
4059          ELSEIF ( lug_palm == ind_luv_de_shrubs  )  THEN
4060             lug_dep = 14
4061          ELSEIF ( lug_palm == ind_luv_mixed_forest )  THEN
4062             lug_dep = 4
4063          ELSEIF ( lug_palm == ind_luv_intrup_forest )  THEN
4064             lug_dep = 8     
4065          ENDIF
4066       ENDIF
4067
4068       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
4069!
4070!--       For walls in USM assume concrete walls/roofs,
4071!--       assumed LU class desert as also assumed for
4072!--       pavements in LSM
4073          luu_palm = ind_lup_conc
4074          IF ( luu_palm == ind_lup_user )  THEN
4075             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
4076             CALL message( 'chem_depo', 'CM0455', 1, 2, 0, 6, 0 )
4077          ELSEIF ( luu_palm == ind_lup_asph_conc )  THEN
4078             luu_dep = 9
4079          ELSEIF ( luu_palm == ind_lup_asph )  THEN
4080             luu_dep = 9
4081          ELSEIF ( luu_palm ==  ind_lup_conc )  THEN
4082             luu_dep = 9
4083          ELSEIF ( luu_palm ==  ind_lup_sett )  THEN
4084             luu_dep = 9
4085          ELSEIF ( luu_palm == ind_lup_pav_stones )  THEN
4086             luu_dep = 9
4087          ELSEIF ( luu_palm == ind_lup_cobblest )  THEN
4088             luu_dep = 9       
4089          ELSEIF ( luu_palm == ind_lup_metal )  THEN
4090             luu_dep = 9
4091          ELSEIF ( luu_palm == ind_lup_wood )  THEN
4092             luu_dep = 9   
4093          ELSEIF ( luu_palm == ind_lup_gravel )  THEN
4094             luu_dep = 9
4095          ELSEIF ( luu_palm == ind_lup_f_gravel )  THEN
4096             luu_dep = 9
4097          ELSEIF ( luu_palm == ind_lup_pebblest )  THEN
4098             luu_dep = 9
4099          ELSEIF ( luu_palm == ind_lup_woodchips )  THEN
4100             luu_dep = 9
4101          ELSEIF ( luu_palm == ind_lup_tartan )  THEN
4102             luu_dep = 9
4103          ELSEIF ( luu_palm == ind_lup_art_turf )  THEN
4104             luu_dep = 9
4105          ELSEIF ( luu_palm == ind_lup_clay )  THEN
4106             luu_dep = 9
4107          ENDIF
4108       ENDIF
4109
4110       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
4111!
4112!--       For windows in USM assume metal as this is
4113!--       as close as we get, assumed LU class desert
4114!--       as also assumed for pavements in LSM
4115          lud_palm = ind_lup_metal     
4116          IF ( lud_palm == ind_lup_user )  THEN
4117             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
4118             CALL message( 'chem_depo', 'CM0456', 1, 2, 0, 6, 0 )
4119          ELSEIF ( lud_palm == ind_lup_asph_conc )  THEN
4120             lud_dep = 9
4121          ELSEIF ( lud_palm == ind_lup_asph )  THEN
4122             lud_dep = 9
4123          ELSEIF ( lud_palm ==  ind_lup_conc )  THEN
4124             lud_dep = 9
4125          ELSEIF ( lud_palm ==  ind_lup_sett )  THEN
4126             lud_dep = 9
4127          ELSEIF ( lud_palm == ind_lup_pav_stones )  THEN
4128             lud_dep = 9
4129          ELSEIF ( lud_palm == ind_lup_cobblest )  THEN
4130             lud_dep = 9       
4131          ELSEIF ( lud_palm == ind_lup_metal )  THEN
4132             lud_dep = 9
4133          ELSEIF ( lud_palm == ind_lup_wood )  THEN
4134             lud_dep = 9   
4135          ELSEIF ( lud_palm == ind_lup_gravel )  THEN
4136             lud_dep = 9
4137          ELSEIF ( lud_palm == ind_lup_f_gravel )  THEN
4138             lud_dep = 9
4139          ELSEIF ( lud_palm == ind_lup_pebblest )  THEN
4140             lud_dep = 9
4141          ELSEIF ( lud_palm == ind_lup_woodchips )  THEN
4142             lud_dep = 9
4143          ELSEIF ( lud_palm == ind_lup_tartan )  THEN
4144             lud_dep = 9
4145          ELSEIF ( lud_palm == ind_lup_art_turf )  THEN
4146             lud_dep = 9
4147          ELSEIF ( lud_palm == ind_lup_clay )  THEN
4148             lud_dep = 9
4149          ENDIF
4150       ENDIF
4151!
4152!--    @TODO: Activate these lines as soon as new ebsolver branch is merged:
4153!--    Set wetness indicator to dry or wet for usm vegetation or pavement
4154       !IF ( surf_usm_h%c_liq(m) > 0 )  THEN
4155       !   nwet = 1
4156       !ELSE
4157       nwet = 0
4158       !ENDIF
4159!
4160!--    Compute length of time step
4161       IF ( call_chem_at_all_substeps )  THEN
4162          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
4163       ELSE
4164          dt_chem = dt_3d
4165       ENDIF
4166
4167       dh = dzw(k)
4168       inv_dh = 1.0_wp / dh
4169       dt_dh = dt_chem/dh
4170!
4171!--    Concentration at i,j,k
4172       DO  lsp = 1, nspec
4173          conc_ijk(lsp) = chem_species(lsp)%conc(k,j,i)
4174       ENDDO
4175!
4176!--    Temperature at i,j,k
4177       temp_tmp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
4178
4179       ts       = temp_tmp - 273.15  !< in degrees celcius
4180!
4181!--    Viscosity of air
4182       visc = 1.496e-6 * temp_tmp**1.5 / (temp_tmp + 120.0)
4183!
4184!--    Air density at k
4185       dens = rho_air_zw(k)
4186!
4187!--    Calculate relative humidity from specific humidity for DEPAC
4188       qv_tmp = MAX(q(k,j,i),0.0_wp)
4189       rh_surf = relativehumidity_from_specifichumidity(qv_tmp, temp_tmp, hyp(k) )
4190!
4191!--    Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
4192!--    for each surface fraction. Then derive overall budget taking into account the surface fractions.
4193!
4194!--    Walls/roofs
4195       IF ( surf_usm_h%frac(ind_veg_wall,m) > 0 )  THEN
4196!
4197!--       No vegetation on non-green walls:
4198          lai = 0.0_wp
4199          sai = 0.0_wp
4200
4201          slinnfac = 1.0_wp
4202!
4203!--       Get vd
4204          DO  lsp = 1, nvar
4205!
4206!--          Initialize
4207             vs = 0.0_wp
4208             vd_lu = 0.0_wp
4209             rs = 0.0_wp
4210             rb = 0.0_wp
4211             rc_tot = 0.0_wp
4212             IF (spc_names(lsp) == 'PM10' )  THEN
4213                part_type = 1
4214!
4215!--             Sedimentation velocity
4216                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4217                     particle_pars(ind_p_size, part_type), &
4218                     particle_pars(ind_p_slip, part_type), &
4219                     visc)
4220
4221                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4222                     vs, &
4223                     particle_pars(ind_p_size, part_type), &
4224                     particle_pars(ind_p_slip, part_type), &
4225                     nwet, temp_tmp, dens, visc, &
4226                     luu_dep,  &
4227                     r_aero_surf, ustar_surf )
4228
4229                bud_luu(lsp) = - conc_ijk(lsp) * &
4230                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4231
4232             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4233                part_type = 2
4234!
4235!--             Sedimentation velocity
4236                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4237                     particle_pars(ind_p_size, part_type), &
4238                     particle_pars(ind_p_slip, part_type), &
4239                     visc)
4240
4241                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4242                     vs, &
4243                     particle_pars(ind_p_size, part_type), &
4244                     particle_pars(ind_p_slip, part_type), &
4245                     nwet, temp_tmp, dens, visc, &
4246                     luu_dep , &
4247                     r_aero_surf, ustar_surf )
4248
4249                bud_luu(lsp) = - conc_ijk(lsp) * &
4250                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4251
4252             ELSE  !< GASES
4253!
4254!--             Read spc_name of current species for gas parameter
4255
4256                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4257                   i_pspec = 0
4258                   DO  pspec = 1, nposp
4259                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4260                         i_pspec = pspec
4261                      END IF
4262                   ENDDO
4263                ELSE
4264!
4265!--                For now species not deposited
4266                   CYCLE
4267                ENDIF
4268!
4269!--             Factor used for conversion from ppb to ug/m3 :
4270!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4271!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4272!--                 c           1e-9              xm_tracer         1e9       /       xm_air            dens
4273!--             thus:
4274!--                 c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4275!--             Use density at k:
4276
4277                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
4278
4279                !
4280!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4281!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4282                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4283!
4284!--             Diffusivity for DEPAC relevant gases
4285!--             Use default value
4286                diffusivity            = 0.11e-4
4287!
4288!--             overwrite with known coefficients of diffusivity from Massman (1998)
4289                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4290                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4291                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4292                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4293                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4294                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4295                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4296!
4297!--             Get quasi-laminar boundary layer resistance rb:
4298                CALL get_rb_cell( (luu_dep == ilu_water_sea) .OR. (luu_dep == ilu_water_inland),   &
4299                     z0h_surf, ustar_surf, diffusivity, &
4300                     rb )
4301!
4302!--             Get rc_tot
4303                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,          &
4304                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, luu_dep, 2,     &
4305                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,   &
4306                                         r_aero_surf, rb )
4307!
4308!--             Calculate budget
4309                IF ( rc_tot <= 0.0 )  THEN
4310
4311                   bud_luu(lsp) = 0.0_wp
4312
4313                ELSE
4314
4315                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4316
4317                   bud_luu(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4318                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4319                ENDIF
4320
4321             ENDIF
4322          ENDDO
4323       ENDIF
4324!
4325!--    Green usm surfaces
4326       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
4327
4328!
4329!--       No vegetation on bare soil, desert or ice:
4330          IF ( ( lug_palm == ind_luv_b_soil ) .OR. &
4331                ( lug_palm == ind_luv_desert ) .OR. &
4332                 ( lug_palm == ind_luv_ice ) ) THEN
4333
4334             lai = 0.0_wp
4335             sai = 0.0_wp
4336
4337          ENDIF
4338
4339         
4340          slinnfac = 1.0_wp
4341!
4342!--       Get vd
4343          DO  lsp = 1, nvar
4344!
4345!--          Initialize
4346             vs = 0.0_wp
4347             vd_lu = 0.0_wp
4348             rs = 0.0_wp
4349             rb = 0.0_wp
4350             rc_tot = 0.0_wp
4351             IF ( spc_names(lsp) == 'PM10' )  THEN
4352                part_type = 1
4353!
4354!--             Sedimentation velocity
4355                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4356                     particle_pars(ind_p_size, part_type), &
4357                     particle_pars(ind_p_slip, part_type), &
4358                     visc)
4359
4360                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4361                     vs, &
4362                     particle_pars(ind_p_size, part_type), &
4363                     particle_pars(ind_p_slip, part_type), &
4364                     nwet, temp_tmp, dens, visc, &
4365                     lug_dep,  &
4366                     r_aero_surf, ustar_surf )
4367
4368                bud_lug(lsp) = - conc_ijk(lsp) * &
4369                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4370
4371             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4372                part_type = 2
4373!
4374!--             Sedimentation velocity
4375                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4376                     particle_pars(ind_p_size, part_type), &
4377                     particle_pars(ind_p_slip, part_type), &
4378                     visc)
4379
4380                CALL drydepo_aero_zhang_vd( vd_lu, rs, &
4381                     vs, &
4382                     particle_pars(ind_p_size, part_type), &
4383                     particle_pars(ind_p_slip, part_type), &
4384                     nwet, temp_tmp, dens, visc, &
4385                     lug_dep, &
4386                     r_aero_surf, ustar_surf )
4387
4388                bud_lug(lsp) = - conc_ijk(lsp) * &
4389                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4390
4391             ELSE  !< GASES
4392!
4393!--             Read spc_name of current species for gas parameter
4394
4395                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4396                   i_pspec = 0
4397                   DO  pspec = 1, nposp
4398                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4399                         i_pspec = pspec
4400                      END IF
4401                   ENDDO
4402                ELSE
4403!
4404!--                For now species not deposited
4405                   CYCLE
4406                ENDIF
4407!
4408!--             Factor used for conversion from ppb to ug/m3 :
4409!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4410!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4411!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
4412!--             thus:
4413!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4414!--             Use density at k:
4415
4416                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
4417!
4418!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4419                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4420                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4421!
4422!--             Diffusivity for DEPAC relevant gases
4423!--             Use default value
4424                diffusivity            = 0.11e-4
4425!
4426!--             overwrite with known coefficients of diffusivity from Massman (1998)
4427                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4428                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4429                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4430                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4431                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4432                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4433                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4434!
4435!--             Get quasi-laminar boundary layer resistance rb:
4436                CALL get_rb_cell( (lug_dep == ilu_water_sea) .OR. (lug_dep == ilu_water_inland),    &
4437                     z0h_surf, ustar_surf, diffusivity, &
4438                     rb )
4439!
4440!--             Get rc_tot
4441                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,           &
4442                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lug_dep, 2,      &
4443                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,    &
4444                                         r_aero_surf , rb )
4445!
4446!--             Calculate budget
4447                IF ( rc_tot <= 0.0 )  THEN
4448
4449                   bud_lug(lsp) = 0.0_wp
4450
4451                ELSE
4452
4453                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4454
4455                   bud_lug(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4456                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4457                ENDIF
4458
4459             ENDIF
4460          ENDDO
4461       ENDIF
4462!
4463!--    Windows
4464       IF ( surf_usm_h%frac(ind_wat_win,m) > 0 )  THEN
4465!
4466!--       No vegetation on windows:
4467          lai = 0.0_wp
4468          sai = 0.0_wp
4469
4470          slinnfac = 1.0_wp
4471!
4472!--       Get vd
4473          DO  lsp = 1, nvar
4474!
4475!--          Initialize
4476             vs = 0.0_wp
4477             vd_lu = 0.0_wp
4478             rs = 0.0_wp
4479             rb = 0.0_wp
4480             rc_tot = 0.0_wp 
4481             IF ( spc_names(lsp) == 'PM10' )  THEN
4482                part_type = 1
4483!
4484!--             Sedimentation velocity
4485                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4486                     particle_pars(ind_p_size, part_type), &
4487                     particle_pars(ind_p_slip, part_type), &
4488                     visc)
4489
4490                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
4491                     particle_pars(ind_p_size, part_type), &
4492                     particle_pars(ind_p_slip, part_type), &
4493                     nwet, temp_tmp, dens, visc,              &
4494                     lud_dep, r_aero_surf, ustar_surf )
4495
4496                bud_lud(lsp) = - conc_ijk(lsp) * &
4497                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4498
4499             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
4500                part_type = 2
4501!
4502!--             Sedimentation velocity
4503                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
4504                     particle_pars(ind_p_size, part_type), &
4505                     particle_pars(ind_p_slip, part_type), &
4506                     visc)
4507
4508                CALL drydepo_aero_zhang_vd( vd_lu, rs, vs, &
4509                     particle_pars(ind_p_size, part_type), &
4510                     particle_pars(ind_p_slip, part_type), &
4511                     nwet, temp_tmp, dens, visc, &
4512                     lud_dep, &
4513                     r_aero_surf, ustar_surf )
4514
4515                bud_lud(lsp) = - conc_ijk(lsp) * &
4516                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4517
4518             ELSE  !< GASES
4519!
4520!--             Read spc_name of current species for gas PARAMETER
4521
4522                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
4523                   i_pspec = 0
4524                   DO  pspec = 1, nposp
4525                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
4526                         i_pspec = pspec
4527                      END IF
4528                   ENDDO
4529                ELSE
4530!
4531!--                For now species not deposited
4532                   CYCLE
4533                ENDIF
4534!
4535!--             Factor used for conversion from ppb to ug/m3 :
4536!--                 ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
4537!--                 (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
4538!--                 c           1e-9               xm_tracer         1e9       /       xm_air            dens
4539!--             thus:
4540!--                  c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
4541!--             Use density at k:
4542
4543                ppm2ugm3 =  (dens/xm_air) * 0.001_wp  ! (mole air)/m3
4544!
4545!--             Atmospheric concentration in DEPAC is requested in ug/m3:
4546!--                 ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
4547                conc_ijk_ugm3 = conc_ijk(lsp)         * ppm2ugm3 *   specmolm(i_pspec)  ! in ug/m3
4548!
4549!--             Diffusivity for DEPAC relevant gases
4550!--             Use default value
4551                diffusivity = 0.11e-4
4552!
4553!--             overwrite with known coefficients of diffusivity from Massman (1998)
4554                IF ( spc_names(lsp) == 'NO2' ) diffusivity = 0.136e-4 
4555                IF ( spc_names(lsp) == 'NO'  ) diffusivity = 0.199e-4
4556                IF ( spc_names(lsp) == 'O3'  ) diffusivity = 0.144e-4
4557                IF ( spc_names(lsp) == 'CO'  ) diffusivity = 0.176e-4
4558                IF ( spc_names(lsp) == 'SO2' ) diffusivity = 0.112e-4
4559                IF ( spc_names(lsp) == 'CH4' ) diffusivity = 0.191e-4
4560                IF ( spc_names(lsp) == 'NH3' ) diffusivity = 0.191e-4
4561!
4562!--             Get quasi-laminar boundary layer resistance rb:
4563                CALL get_rb_cell( (lud_dep == ilu_water_sea) .OR. (lud_dep == ilu_water_inland),   &
4564                     z0h_surf, ustar_surf, diffusivity, rb )
4565!
4566!--             Get rc_tot
4567                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_surf,         &
4568                                         solar_rad, cos_zenith, rh_surf, lai, sai, nwet, lud_dep, 2,    &
4569                                         rc_tot, ccomp_tot(lsp), hyp(nzb), conc_ijk_ugm3, diffusivity,  &
4570                                         r_aero_surf , rb )
4571!
4572!--             Calculate budget
4573                IF ( rc_tot <= 0.0 )  THEN
4574
4575                   bud_lud(lsp) = 0.0_wp
4576
4577                ELSE
4578
4579                   vd_lu = 1.0_wp / (r_aero_surf + rb + rc_tot )
4580
4581                   bud_lud(lsp) = - (conc_ijk(lsp) - ccomp_tot(lsp)) * &
4582                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
4583                ENDIF
4584
4585             ENDIF
4586          ENDDO
4587       ENDIF
4588
4589
4590       bud = 0.0_wp
4591!
4592!--    Calculate overall budget for surface m and adapt concentration
4593       DO  lsp = 1, nspec
4594
4595
4596          bud(lsp) = surf_usm_h%frac(ind_veg_wall,m) * bud_luu(lsp) + &
4597               surf_usm_h%frac(ind_pav_green,m) * bud_lug(lsp) + &
4598               surf_usm_h%frac(ind_wat_win,m) * bud_lud(lsp)
4599!
4600!--       Compute new concentration
4601          conc_ijk(lsp) = conc_ijk(lsp) + bud(lsp) * inv_dh
4602
4603          chem_species(lsp)%conc(k,j,i) = MAX( 0.0_wp, conc_ijk(lsp) )
4604
4605       ENDDO
4606
4607    ENDIF
4608
4609
4610 END SUBROUTINE chem_depo
4611
4612
4613!------------------------------------------------------------------------------!
4614! Description:
4615! ------------
4616!> Subroutine to compute total canopy (or surface) resistance Rc for gases
4617!>
4618!> DEPAC:
4619!> Code of the DEPAC routine and corresponding subroutines below from the DEPAC
4620!> module of the LOTOS-EUROS model (Manders et al., 2017)
4621!>
4622!> Original DEPAC routines by RIVM and TNO (2015), for Documentation see
4623!> van Zanten et al., 2010.
4624!------------------------------------------------------------------------------!
4625 SUBROUTINE drydepos_gas_depac( compnam, day_of_year, lat, t, ust, solar_rad, sinphi,    &
4626      rh, lai, sai, nwet, lu, iratns, rc_tot, ccomp_tot, p, conc_ijk_ugm3, diffusivity,  &
4627      ra, rb ) 
4628!
4629!--   Some of depac arguments are OPTIONAL:
4630!--    A. compute Rc_tot without compensation points (ccomp_tot will be zero):
4631!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi])
4632!--    B. compute Rc_tot with compensation points (used for LOTOS-EUROS):
4633!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4634!--                c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water)
4635!--
4636!--    C. compute effective Rc based on compensation points (used for OPS):
4637!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4638!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4639!--               ra, rb, rc_eff)
4640!--    X1. Extra (OPTIONAL) output variables:
4641!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4642!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4643!--               ra, rb, rc_eff, &
4644!--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out)
4645!--    X2. Extra (OPTIONAL) needed for stomatal ozone flux calculation (only sunlit leaves):
4646!--        CALL depac (compnam, day_of_year, lat, t, ust, glrad, sinphi, rh, nwet, lu, iratns, rc_tot, ccomp_tot, [smi], &
4647!--               c_ave_prev_nh3, c_ave_prev_so2, catm, gamma_soil_water, &
4648!--               ra, rb, rc_eff, &
4649!--               gw_out, gstom_out, gsoil_eff_out, cw_out, cstom_out, csoil_out, lai_out, sai_out, &
4650!--               calc_stom_o3flux, frac_sto_o3_lu, fac_surface_area_2_PLA)
4651
4652
4653    CHARACTER(LEN=*), INTENT(IN) ::  compnam         !< component name
4654                                                     !< 'HNO3','NO','NO2','O3','SO2','NH3'
4655    INTEGER(iwp), INTENT(IN) ::  day_of_year         !< day of year, 1 ... 365 (366)
4656    INTEGER(iwp), INTENT(IN) ::  nwet                !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4657    INTEGER(iwp), INTENT(IN) ::  lu                  !< land use type, lu = 1,...,nlu
4658    INTEGER(iwp), INTENT(IN) ::  iratns              !< index for NH3/SO2 ratio used for SO2:
4659                                                     !< iratns = 1: low NH3/SO2
4660                                                     !< iratns = 2: high NH3/SO2
4661                                                     !< iratns = 3: very low NH3/SO2
4662    REAL(wp), INTENT(IN) ::  lat                     !< latitude Northern hemisphere (degrees) (S. hemisphere not possible)
4663    REAL(wp), INTENT(IN) ::  t                       !< temperature (C)
4664    REAL(wp), INTENT(IN) ::  ust                     !< friction velocity (m/s)
4665    REAL(wp), INTENT(IN) ::  solar_rad               !< solar radiation, dirict+diffuse (W/m2)
4666    REAL(wp), INTENT(IN) ::  sinphi                  !< sin of solar elevation angle
4667    REAL(wp), INTENT(IN) ::  rh                      !< relative humidity (%)
4668    REAL(wp), INTENT(IN) ::  lai                     !< one-sidedleaf area index (-)
4669    REAL(wp), INTENT(IN) ::  sai                     !< surface area index (-) (lai + branches and stems)
4670    REAL(wp), INTENT(IN) ::  diffusivity             !< diffusivity
4671    REAL(wp), INTENT(IN) ::  p                       !< pressure (Pa)
4672    REAL(wp), INTENT(IN) ::  conc_ijk_ugm3           !< actual atmospheric concentration (ug/m3), in DEPAC=Catm
4673    REAL(wp), INTENT(IN) ::  ra                      !< aerodynamic resistance (s/m)
4674    REAL(wp), INTENT(IN) ::  rb                      !< boundary layer resistance (s/m)
4675
4676    REAL(wp), INTENT(OUT) ::  rc_tot                 !< total canopy resistance Rc (s/m)
4677    REAL(wp), INTENT(OUT) ::  ccomp_tot              !< total compensation point (ug/m3)
4678!                                                     !< [= 0 for species that don't have a compensation
4679!-- Local variables:
4680!
4681!-- Component number taken from component name, paramteres matched with include files
4682    INTEGER(iwp) ::  icmp
4683!
4684!-- Component numbers:
4685    INTEGER(iwp), PARAMETER ::  icmp_o3   = 1
4686    INTEGER(iwp), PARAMETER ::  icmp_so2  = 2
4687    INTEGER(iwp), PARAMETER ::  icmp_no2  = 3
4688    INTEGER(iwp), PARAMETER ::  icmp_no   = 4
4689    INTEGER(iwp), PARAMETER ::  icmp_nh3  = 5
4690    INTEGER(iwp), PARAMETER ::  icmp_co   = 6
4691    INTEGER(iwp), PARAMETER ::  icmp_no3  = 7
4692    INTEGER(iwp), PARAMETER ::  icmp_hno3 = 8
4693    INTEGER(iwp), PARAMETER ::  icmp_n2o5 = 9
4694    INTEGER(iwp), PARAMETER ::  icmp_h2o2 = 10
4695
4696    LOGICAL ::  ready                                !< Rc has been set:
4697    !< = 1 -> constant Rc
4698    !< = 2 -> temperature dependent Rc
4699!
4700!-- Vegetation indicators:
4701    LOGICAL ::  lai_present                          !< leaves are present for current land use type
4702    LOGICAL ::  sai_present                          !< vegetation is present for current land use type
4703
4704!    REAL(wp) ::  laimax                              !< maximum leaf area index (-)
4705    REAL(wp) ::  gw                                  !< external leaf conductance (m/s)
4706    REAL(wp) ::  gstom                               !< stomatal conductance (m/s)
4707    REAL(wp) ::  gsoil_eff                           !< effective soil conductance (m/s)
4708    REAL(wp) ::  gc_tot                              !< total canopy conductance (m/s)
4709    REAL(wp) ::  cw                                  !< external leaf surface compensation point (ug/m3)
4710    REAL(wp) ::  cstom                               !< stomatal compensation point (ug/m3)
4711    REAL(wp) ::  csoil                               !< soil compensation point (ug/m3)
4712!
4713!-- Next statement is just to avoid compiler warning about unused variable
4714    IF ( day_of_year == 0  .OR.  ( conc_ijk_ugm3 + lat + ra + rb ) > 0.0_wp )  CONTINUE
4715!
4716!-- Define component number
4717    SELECT CASE ( TRIM( compnam ) )
4718
4719    CASE ( 'O3', 'o3' )
4720       icmp = icmp_o3
4721
4722    CASE ( 'SO2', 'so2' )
4723       icmp = icmp_so2
4724
4725    CASE ( 'NO2', 'no2' )
4726       icmp = icmp_no2
4727
4728    CASE ( 'NO', 'no' )
4729       icmp = icmp_no
4730
4731    CASE ( 'NH3', 'nh3' )
4732       icmp = icmp_nh3
4733
4734    CASE ( 'CO', 'co' )
4735       icmp = icmp_co
4736
4737    CASE ( 'NO3', 'no3' )
4738       icmp = icmp_no3
4739
4740    CASE ( 'HNO3', 'hno3' )
4741       icmp = icmp_hno3
4742
4743    CASE ( 'N2O5', 'n2o5' )
4744       icmp = icmp_n2o5
4745
4746    CASE ( 'H2O2', 'h2o2' )
4747       icmp = icmp_h2o2
4748
4749    CASE default
4750!
4751!--    Component not part of DEPAC --> not deposited
4752       RETURN
4753
4754    END SELECT
4755
4756!
4757!-- Inititalize
4758    gw        = 0.0_wp
4759    gstom     = 0.0_wp
4760    gsoil_eff = 0.0_wp
4761    gc_tot    = 0.0_wp
4762    cw        = 0.0_wp
4763    cstom     = 0.0_wp
4764    csoil     = 0.0_wp
4765!
4766!-- Check whether vegetation is present:
4767    lai_present = ( lai > 0.0 )
4768    sai_present = ( sai > 0.0 )
4769!
4770!-- Set Rc (i.e. rc_tot) in special cases:
4771    CALL rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
4772!
4773!-- If Rc is not set:
4774    IF ( .NOT. ready ) then
4775!
4776!--    External conductance:
4777       CALL rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai,gw )         
4778!
4779!--    Stomatal conductance:
4780       CALL rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p )
4781!
4782!--    Effective soil conductance:
4783       CALL rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff )
4784!
4785!--    Total canopy conductance (gc_tot) and resistance Rc (rc_tot):
4786       CALL rc_rctot( gstom, gsoil_eff, gw, gc_tot, rc_tot )
4787!
4788!--    Needed to include compensation point for NH3
4789!--    Compensation points (always returns ccomp_tot; currently ccomp_tot != 0 only for NH3):
4790!--    CALL rc_comp_point( compnam,lu,day_of_year,t,gw,gstom,gsoil_eff,gc_tot,&
4791!--          lai_present, sai_present, &
4792!--          ccomp_tot, &
4793!--          conc_ijk_ugm3=conc_ijk_ugm3,c_ave_prev_nh3=c_ave_prev_nh3, &
4794!--          c_ave_prev_so2=c_ave_prev_so2,gamma_soil_water=gamma_soil_water, &
4795!--          tsea=tsea,cw=cw,cstom=cstom,csoil=csoil )
4796!
4797!--    Effective Rc based on compensation points:
4798!--        IF ( present(rc_eff) ) then
4799!--          check on required arguments:
4800!--           IF ( (.not. present(conc_ijk_ugm3)) .OR. (.not. present(ra)) .OR. (.not. present(rb)) ) then
4801!--              stop 'output argument rc_eff requires input arguments conc_ijk_ugm3, ra and rb'
4802!--           END IF
4803!
4804!--       Compute rc_eff :
4805       !      CALL rc_comp_point_rc_eff(ccomp_tot,conc_ijk_ugm3,ra,rb,rc_tot,rc_eff)
4806       !   ENDIF
4807    ENDIF
4808
4809 END SUBROUTINE drydepos_gas_depac
4810
4811
4812!------------------------------------------------------------------------------!
4813! Description:
4814! ------------
4815!> Subroutine to compute total canopy resistance in special cases
4816!------------------------------------------------------------------------------!
4817 SUBROUTINE rc_special( icmp, compnam, lu, t, nwet, rc_tot, ready, ccomp_tot )
4818
4819   
4820    CHARACTER(LEN=*), INTENT(IN)  ::  compnam     !< component name
4821
4822    INTEGER(iwp), INTENT(IN)  ::  icmp            !< component index     
4823    INTEGER(iwp), INTENT(IN)  ::  lu              !< land use type, lu = 1,...,nlu
4824    INTEGER(iwp), INTENT(IN)  ::  nwet            !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4825
4826    REAL(wp), INTENT(IN)  ::  t                   !< temperature (C)
4827
4828    REAL(wp), INTENT(OUT) ::  rc_tot             !< total canopy resistance Rc (s/m)
4829    REAL(wp), INTENT(OUT) ::  ccomp_tot          !< total compensation point (ug/m3)
4830
4831    LOGICAL, INTENT(OUT) ::  ready               !< Rc has been set
4832                                                 !< = 1 -> constant Rc
4833!
4834!-- Next line is to avoid compiler warning about unused variable
4835    IF ( icmp == 0 )  CONTINUE
4836!
4837!-- rc_tot is not yet set:
4838    ready = .FALSE.
4839!
4840!-- Default compensation point in special CASEs = 0:
4841    ccomp_tot = 0.0_wp
4842
4843    SELECT CASE( TRIM( compnam ) )
4844    CASE( 'HNO3', 'N2O5', 'NO3', 'H2O2' )
4845!
4846!--    No separate resistances for HNO3; just one total canopy resistance:
4847       IF ( t < -5.0_wp .AND. nwet == 9 )  THEN
4848!
4849!--       T < 5 C and snow:
4850          rc_tot = 50.0_wp
4851       ELSE
4852!
4853!--       all other circumstances:
4854          rc_tot = 10.0_wp
4855       ENDIF
4856       ready = .TRUE.
4857
4858    CASE( 'NO', 'CO' )
4859       IF ( lu == ilu_water_sea .OR. lu == ilu_water_inland )  THEN       ! water
4860          rc_tot = 2000.0_wp
4861          ready = .TRUE.
4862       ELSEIF ( nwet == 1 )  THEN  !< wet
4863          rc_tot = 2000.0_wp
4864          ready = .TRUE.
4865       ENDIF
4866    CASE( 'NO2', 'O3', 'SO2', 'NH3' )
4867!
4868!--    snow surface:
4869       IF ( nwet == 9 )  THEN
4870!
4871!--       To be activated when snow is implemented
4872          !CALL rc_snow(ipar_snow(icmp),t,rc_tot)
4873          ready = .TRUE.
4874       ENDIF
4875    CASE default
4876       message_string = 'Component '// TRIM( compnam ) // ' not supported'
4877       CALL message( 'rc_special', 'CM0457', 1, 2, 0, 6, 0 )
4878    END SELECT
4879
4880 END SUBROUTINE rc_special
4881
4882
4883!------------------------------------------------------------------------------!
4884! Description:
4885! ------------
4886!> Subroutine to compute external conductance
4887!------------------------------------------------------------------------------!
4888 SUBROUTINE rc_gw( compnam, iratns, t, rh, nwet, sai_present, sai, gw )
4889
4890!
4891!-- Input/output variables:
4892    CHARACTER(LEN=*), INTENT(IN) ::  compnam      !< component name
4893
4894    INTEGER(iwp), INTENT(IN) ::  nwet             !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4895    INTEGER(iwp), INTENT(IN) ::  iratns           !< index for NH3/SO2 ratio;
4896                                                  !< iratns = 1: low NH3/SO2
4897                                                  !< iratns = 2: high NH3/SO2
4898                                                  !< iratns = 3: very low NH3/SO2
4899    LOGICAL, INTENT(IN) ::  sai_present
4900
4901    REAL(wp), INTENT(IN) ::  t                    !< temperature (C)
4902    REAL(wp), INTENT(IN) ::  rh                   !< relative humidity (%)
4903    REAL(wp), INTENT(IN) ::  sai                  !< one-sided leaf area index (-)
4904
4905    REAL(wp), INTENT(OUT) ::  gw                  !< external leaf conductance (m/s)
4906
4907    SELECT CASE( TRIM( compnam ) )
4908
4909    CASE( 'NO2' )
4910       CALL rw_constant( 2000.0_wp, sai_present, gw )
4911
4912    CASE( 'NO', 'CO' )
4913       CALL rw_constant( -9999.0_wp, sai_present, gw )   !< see Erisman et al, 1994 section 3.2.3
4914
4915    CASE( 'O3' )
4916       CALL rw_constant( 2500.0_wp, sai_present, gw )
4917
4918    CASE( 'SO2' )
4919       CALL rw_so2( t, nwet, rh, iratns, sai_present, gw )
4920
4921    CASE( 'NH3' )
4922       CALL rw_nh3_sutton( t, rh, sai_present, gw )
4923!
4924!--    conversion from leaf resistance to canopy resistance by multiplying with sai:
4925       gw = sai * gw
4926
4927    CASE default
4928       message_string = 'Component '// TRIM( compnam ) // ' not supported'
4929       CALL message( 'rc_gw', 'CM0458', 1, 2, 0, 6, 0 )
4930    END SELECT
4931
4932 END SUBROUTINE rc_gw
4933
4934
4935!------------------------------------------------------------------------------!
4936! Description:
4937! ------------
4938!> Subroutine to compute external leaf conductance for SO2
4939!------------------------------------------------------------------------------!
4940 SUBROUTINE rw_so2( t, nwet, rh, iratns, sai_present, gw )
4941
4942!
4943!-- Input/output variables:
4944    INTEGER(iwp), INTENT(IN) ::  nwet        !< wetness indicator; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
4945    INTEGER(iwp), INTENT(IN) ::  iratns      !< index for NH3/SO2 ratio:
4946                                             !< iratns = 1: low NH3/SO2
4947                                             !< iratns = 2: high NH3/SO2
4948                                             !< iratns = 3: very low NH3/SO2
4949    LOGICAL, INTENT(IN) ::  sai_present
4950
4951    REAL(wp), INTENT(IN) ::  t               !< temperature (C)
4952    REAL(wp), INTENT(IN) ::  rh              !< relative humidity (%)   
4953
4954    REAL(wp), INTENT(OUT) ::  gw             !< external leaf conductance (m/s)
4955!
4956!-- Local variables:
4957    REAL(wp) ::  rw                          !< external leaf resistance (s/m)
4958!
4959!-- Check if vegetation present:
4960    IF ( sai_present )  THEN
4961
4962       IF ( nwet == 0 )  THEN
4963!
4964!--   ------------------------
4965!--         dry surface
4966!--   ------------------------
4967!--         T > -1 C
4968          IF ( t > -1.0_wp )  THEN
4969             IF ( rh < 81.3_wp )  THEN
4970                rw = 25000.0_wp * exp( -0.0693_wp * rh )
4971             ELSE
4972                rw = 0.58e12 * exp( -0.278_wp * rh ) + 10.0_wp
4973             ENDIF
4974          ELSE
4975             ! -5 C < T <= -1 C
4976             IF ( t > -5.0_wp )  THEN
4977                rw = 200.0_wp
4978             ELSE
4979                ! T <= -5 C
4980                rw = 500.0_wp
4981             ENDIF
4982          ENDIF
4983       ELSE
4984!
4985!--   ------------------------
4986!--         wet surface
4987!--   ------------------------
4988          rw = 10.0_wp !see Table 5, Erisman et al, 1994 Atm. Environment, 0 is impl. as 10
4989       ENDIF
4990!
4991!--    very low NH3/SO2 ratio:
4992       IF ( iratns == iratns_very_low ) rw = rw + 50.0_wp
4993!
4994!--      Conductance:
4995       gw = 1.0_wp / rw
4996    ELSE
4997!
4998!--      no vegetation:
4999       gw = 0.0_wp
5000    ENDIF
5001
5002 END SUBROUTINE rw_so2
5003
5004
5005!------------------------------------------------------------------------------!
5006! Description:
5007! ------------
5008!> Subroutine to compute external leaf conductance for NH3,
5009!>                  following Sutton & Fowler, 1993
5010!------------------------------------------------------------------------------!
5011 SUBROUTINE rw_nh3_sutton( tsurf, rh,sai_present, gw )
5012
5013!
5014!-- Input/output variables:
5015    LOGICAL, INTENT(IN) ::  sai_present
5016
5017    REAL(wp), INTENT(IN) ::  tsurf          !< surface temperature (C)
5018    REAL(wp), INTENT(IN) ::  rh             !< relative humidity (%)
5019
5020    REAL(wp), INTENT(OUT) ::  gw            !< external leaf conductance (m/s)
5021!
5022!-- Local variables:
5023    REAL(wp) ::  rw                         !< external leaf resistance (s/m)
5024    REAL(wp) ::  sai_grass_haarweg          !< surface area index at experimental site Haarweg
5025!
5026!-- Fix sai_grass at value valid for Haarweg data for which gamma_w parametrization is derived
5027    sai_grass_haarweg = 3.5_wp
5028!
5029!-- Calculation rw:
5030!--                    100 - rh
5031!--    rw = 2.0 * exp(----------)
5032!--                      12
5033
5034    IF ( sai_present )  THEN
5035!
5036!--    External resistance according to Sutton & Fowler, 1993
5037       rw = 2.0_wp * exp( ( 100.0_wp - rh ) / 12.0_wp )
5038       rw = sai_grass_haarweg * rw
5039!
5040!--    Frozen soil (from Depac v1):
5041       IF ( tsurf < 0.0_wp ) rw = 200.0_wp
5042!
5043!--    Conductance:
5044       gw = 1.0_wp / rw
5045    ELSE
5046       ! no vegetation:
5047       gw = 0.0_wp
5048    ENDIF
5049
5050 END SUBROUTINE rw_nh3_sutton
5051
5052
5053!------------------------------------------------------------------------------!
5054! Description:
5055! ------------
5056!> Subroutine to compute external leaf conductance
5057!------------------------------------------------------------------------------!
5058 SUBROUTINE rw_constant( rw_val, sai_present, gw )
5059
5060!
5061!-- Input/output variables:
5062    LOGICAL, INTENT(IN) ::  sai_present
5063
5064    REAL(wp), INTENT(IN) ::  rw_val       !< constant value of Rw   
5065
5066    REAL(wp), INTENT(OUT) ::  gw          !< wernal leaf conductance (m/s)
5067!
5068!-- Compute conductance:
5069    IF ( sai_present .AND. .NOT.missing(rw_val) )  THEN
5070       gw = 1.0_wp / rw_val
5071    ELSE
5072       gw = 0.0_wp
5073    ENDIF
5074
5075 END SUBROUTINE rw_constant
5076
5077
5078!------------------------------------------------------------------------------!
5079! Description:
5080! ------------
5081!> Subroutine to compute stomatal conductance
5082!------------------------------------------------------------------------------!
5083 SUBROUTINE rc_gstom( icmp, compnam, lu, lai_present, lai, solar_rad, sinphi, t, rh, diffusivity, gstom, p ) 
5084
5085!
5086!-- input/output variables:
5087    CHARACTER(LEN=*), INTENT(IN) ::  compnam       !< component name
5088
5089    INTEGER(iwp), INTENT(IN) ::  icmp              !< component index
5090    INTEGER(iwp), INTENT(IN) ::  lu                !< land use type , lu = 1,...,nlu
5091
5092    LOGICAL, INTENT(IN) ::  lai_present
5093
5094    REAL(wp), INTENT(IN) ::  lai                   !< one-sided leaf area index
5095    REAL(wp), INTENT(IN) ::  solar_rad             !< solar radiation, dirict+diffuse (W/m2)
5096    REAL(wp), INTENT(IN) ::  sinphi                !< sin of solar elevation angle
5097    REAL(wp), INTENT(IN) ::  t                     !< temperature (C)
5098    REAL(wp), INTENT(IN) ::  rh                    !< relative humidity (%)
5099    REAL(wp), INTENT(IN) ::  diffusivity           !< diffusion coefficient of the gas involved
5100
5101    REAL(wp), OPTIONAL,INTENT(IN) :: p             !< pressure (Pa)
5102
5103    REAL(wp), INTENT(OUT) ::  gstom                !< stomatal conductance (m/s)
5104!
5105!-- Local variables
5106    REAL(wp) ::  vpd                               !< vapour pressure deficit (kPa)
5107
5108    REAL(wp), PARAMETER ::  dO3 = 0.13e-4          !< diffusion coefficient of ozon (m2/s)
5109!
5110!-- Next line is to avoid compiler warning about unused variables
5111    IF ( icmp == 0 )  CONTINUE
5112
5113    SELECT CASE( TRIM( compnam ) )
5114
5115    CASE( 'NO', 'CO' )
5116!
5117!--    For no stomatal uptake is neglected:
5118       gstom = 0.0_wp
5119
5120    CASE( 'NO2', 'O3', 'SO2', 'NH3' )
5121!
5122!--    if vegetation present:
5123       IF ( lai_present )  THEN
5124
5125          IF ( solar_rad > 0.0_wp )  THEN
5126             CALL rc_get_vpd( t, rh, vpd )
5127             CALL rc_gstom_emb( lu, solar_rad, t, vpd, lai_present, lai, sinphi, gstom, p )
5128             gstom = gstom * diffusivity / dO3       !< Gstom of Emberson is derived for ozone
5129          ELSE
5130             gstom = 0.0_wp
5131          ENDIF
5132       ELSE
5133!
5134!--       no vegetation; zero conductance (infinite resistance):
5135          gstom = 0.0_wp
5136       ENDIF
5137
5138    CASE default
5139       message_string = 'Component '// TRIM( compnam ) // ' not supported'
5140       CALL message( 'rc_gstom', 'CM0459', 1, 2, 0, 6, 0 )
5141    END SELECT
5142
5143 END SUBROUTINE rc_gstom
5144
5145
5146!------------------------------------------------------------------------------!
5147! Description:
5148! ------------
5149!> Subroutine to compute stomatal conductance according to Emberson
5150!------------------------------------------------------------------------------!
5151 SUBROUTINE rc_gstom_emb( lu, solar_rad, T, vpd, lai_present, lai, sinp, Gsto, p )
5152!
5153!>  History
5154!>   Original code from Lotos-Euros, TNO, M. Schaap
5155!>   2009-08, M.C. van Zanten, Rivm
5156!>     Updated and extended.
5157!>   2009-09, Arjo Segers, TNO
5158!>     Limitted temperature influence to range to avoid
5159!>     floating point exceptions.
5160
5161!> Method
5162
5163!>   Code based on Emberson et al, 2000, Env. Poll., 403-413
5164!>   Notation conform Unified EMEP Model Description Part 1, ch 8
5165!
5166!>   In the calculation of f_light the modification of L. Zhang 2001, AE to the PARshade and PARsun
5167!>   parametrizations of Norman 1982 are applied
5168!>   f_phen and f_SWP are set to 1
5169!
5170!>   Land use types DEPAC versus Emberson (Table 5.1, EMEP model description)
5171!>   DEPAC                     Emberson
5172!>     1 = grass                 GR = grassland
5173!>     2 = arable land           TC = temperate crops ( lai according to RC = rootcrops)
5174!>     3 = permanent crops       TC = temperate crops ( lai according to RC = rootcrops)
5175!>     4 = coniferous forest     CF = tempareate/boREAL(wp) coniferous forest
5176!>     5 = deciduous forest      DF = temperate/boREAL(wp) deciduous forest
5177!>     6 = water                 W  = water
5178!>     7 = urban                 U  = urban
5179!>     8 = other                 GR = grassland
5180!>     9 = desert                DE = desert
5181!
5182!-- Emberson specific declarations
5183!
5184!-- Input/output variables:
5185    INTEGER(iwp), INTENT(IN) ::  lu             !< land use type, lu = 1,...,nlu
5186
5187    LOGICAL, INTENT(IN) ::  lai_present
5188
5189    REAL(wp), INTENT(IN) ::  solar_rad          !< solar radiation, dirict+diffuse (W/m2)
5190    REAL(wp), INTENT(IN) ::  t                  !< temperature (C)
5191    REAL(wp), INTENT(IN) ::  vpd                !< vapour pressure deficit (kPa)
5192
5193    REAL(wp), INTENT(IN) ::  lai                !< one-sided leaf area index
5194    REAL(wp), INTENT(IN) ::  sinp               !< sin of solar elevation angle
5195
5196    REAL(wp), OPTIONAL, INTENT(IN) ::  p        !< pressure (Pa)
5197
5198    REAL(wp), INTENT(OUT) ::  gsto              !< stomatal conductance (m/s)
5199!
5200!-- Local variables:
5201    REAL(wp) ::  f_light
5202    REAL(wp) ::  f_phen
5203    REAL(wp) ::  f_temp
5204    REAL(wp) ::  f_vpd
5205    REAL(wp) ::  f_swp
5206    REAL(wp) ::  bt
5207    REAL(wp) ::  f_env
5208    REAL(wp) ::  pardir
5209    REAL(wp) ::  pardiff
5210    REAL(wp) ::  parshade
5211    REAL(wp) ::  parsun
5212    REAL(wp) ::  laisun
5213    REAL(wp) ::  laishade
5214    REAL(wp) ::  sinphi
5215    REAL(wp) ::  pres
5216    REAL(wp), PARAMETER ::  p_sealevel = 1.01325e05    !< Pa
5217!
5218!-- Check whether vegetation is present:
5219    IF ( lai_present )  THEN
5220
5221       ! calculation of correction factors for stomatal conductance
5222       IF ( sinp <= 0.0_wp )  THEN 
5223          sinphi = 0.0001_wp
5224       ELSE
5225          sinphi = sinp
5226       END IF
5227!
5228!--    ratio between actual and sea-level pressure is used
5229!--    to correct for height in the computation of par;
5230!--    should not exceed sea-level pressure therefore ...
5231       IF (  present(p) )  THEN
5232          pres = min( p, p_sealevel )
5233       ELSE
5234          pres = p_sealevel
5235       ENDIF
5236!
5237!--    direct and diffuse par, Photoactive (=visible) radiation:
5238       CALL par_dir_diff( solar_rad, sinphi, pres, p_sealevel, pardir, pardiff )
5239!
5240!--    par for shaded leaves (canopy averaged):
5241       parshade = pardiff * exp( -0.5 * lai**0.7 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )     !< Norman,1982
5242       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
5243          parshade = pardiff * exp( -0.5 * lai**0.8 ) + 0.07 * pardir * ( 1.1 - 0.1 * lai ) * exp( -sinphi )  !< Zhang et al., 2001
5244       END IF
5245!
5246!--    par for sunlit leaves (canopy averaged):
5247!--    alpha -> mean angle between leaves and the sun is fixed at 60 deg -> i.e. cos alpha = 0.5
5248       parsun = pardir * 0.5/sinphi + parshade             !< Norman, 1982
5249       IF ( solar_rad > 200.0_wp .AND. lai > 2.5_wp )  THEN
5250          parsun = pardir**0.8 * 0.5 / sinphi + parshade   !< Zhang et al., 2001
5251       END IF
5252!
5253!--    leaf area index for sunlit and shaded leaves:
5254       IF ( sinphi > 0 )  THEN
5255          laisun = 2 * sinphi * ( 1 - exp( -0.5 * lai / sinphi ) )
5256          laishade = lai - laisun
5257       ELSE
5258          laisun = 0
5259          laishade = lai
5260       END IF
5261
5262       f_light = ( laisun * ( 1 - exp( -1.0_wp * alpha(lu) * parsun ) ) + &
5263            laishade * ( 1 - exp( -1.0_wp * alpha(lu) * parshade ) ) ) / lai
5264
5265       f_light = MAX(f_light,f_min(lu))
5266!
5267!--    temperature influence; only non-zero within range [tmin,tmax]:
5268       IF ( ( tmin(lu) < t ) .AND. ( t < tmax(lu) ) )  THEN
5269          bt = ( tmax(lu) - topt(lu) ) / ( topt(lu) - tmin(lu) )
5270          f_temp = ( ( t - tmin(lu) ) / ( topt(lu) - tmin(lu) ) ) * ( ( tmax(lu) - t ) / ( tmax(lu) - topt(lu) ) )**bt
5271       ELSE
5272          f_temp = 0.0_wp
5273       END IF
5274       f_temp = MAX( f_temp, f_min(lu) )
5275!
5276!--      vapour pressure deficit influence
5277       f_vpd = MIN( 1.0_wp, ( ( 1.0_wp - f_min(lu) ) * ( vpd_min(lu) - vpd ) / ( vpd_min(lu) - vpd_max(lu) ) + f_min(lu) ) )
5278       f_vpd = MAX( f_vpd, f_min(lu) )
5279
5280       f_swp = 1.0_wp
5281!
5282!--      influence of phenology on stom. conductance
5283!--      ignored for now in DEPAC since influence of f_phen on lu classes in use is negligible.
5284!--      When other EMEP classes (e.g. med. broadleaf) are used f_phen might be too important to ignore
5285       f_phen = 1.0_wp
5286!
5287!--      evaluate total stomatal conductance
5288       f_env = f_temp * f_vpd * f_swp
5289       f_env = MAX( f_env,f_min(lu) )
5290       gsto = g_max(lu) * f_light * f_phen * f_env
5291!
5292!--      gstom expressed per m2 leafarea;
5293!--      this is converted with lai to m2 surface.
5294       gsto = lai * gsto    ! in m/s
5295
5296    ELSE
5297       gsto = 0.0_wp
5298    ENDIF
5299
5300 END SUBROUTINE rc_gstom_emb
5301
5302
5303 !-------------------------------------------------------------------
5304 !> par_dir_diff
5305 !>     Weiss, A., Norman, J.M. (1985) Partitioning solar radiation into direct and
5306 !>     diffuse, visible and near-infrared components. Agric. Forest Meteorol.
5307 !>     34, 205-213.
5308 !>     From a SUBROUTINE obtained from Leiming Zhang,
5309 !>     Meteorological Service of Canada
5310 !>     Leiming uses solar irradiance. This should be equal to global radiation and
5311 !>     Willem Asman set it to global radiation (here defined as solar radiation, dirict+diffuse)
5312 !>
5313 !>     @todo Check/connect/replace with radiation_model_mod variables   
5314 !-------------------------------------------------------------------
5315 SUBROUTINE par_dir_diff( solar_rad, sinphi, pres, pres_0, par_dir, par_diff )
5316
5317
5318    REAL(wp), INTENT(IN) ::  solar_rad       !< solar radiation, dirict+diffuse (W m-2)
5319    REAL(wp), INTENT(IN) ::  sinphi          !< sine of the solar elevation
5320    REAL(wp), INTENT(IN) ::  pres            !< actual pressure (to correct for height) (Pa)
5321    REAL(wp), INTENT(IN) ::  pres_0          !< pressure at sea level (Pa)
5322
5323    REAL(wp), INTENT(OUT) ::  par_dir        !< par direct : visible (photoactive) direct beam radiation (W m-2)
5324    REAL(wp), INTENT(OUT) ::  par_diff       !< par diffuse: visible (photoactive) diffuse radiation (W m-2)
5325
5326
5327    REAL(wp) ::  sv                          !< total visible radiation
5328    REAL(wp) ::  fv                          !< par direct beam fraction (dimensionless)
5329    REAL(wp) ::  ratio                       !< ratio measured to potential solar radiation (dimensionless)
5330    REAL(wp) ::  rdm                         !< potential direct beam near-infrared radiation (W m-2); "potential" means clear-sky
5331    REAL(wp) ::  rdn                         !< potential diffuse near-infrared radiation (W m-2)
5332    REAL(wp) ::  rdu                         !< visible (par) direct beam radiation (W m-2)
5333    REAL(wp) ::  rdv                         !< potential visible (par) diffuse radiation (W m-2)
5334    REAL(wp) ::  rn                          !< near-infrared radiation (W m-2)
5335    REAL(wp) ::  rv                          !< visible radiation (W m-2)
5336    REAL(wp) ::  ww                          !< water absorption in the near infrared for 10 mm of precipitable water
5337
5338!
5339!-- Calculate visible (PAR) direct beam radiation
5340!-- 600 W m-2 represents average amount of par (400-700 nm wavelength)
5341!-- at the top of the atmosphere; this is roughly 0.45*solar constant (solar constant=1320 Wm-2)
5342    rdu = 600.0_wp* exp( -0.185_wp * ( pres / pres_0 ) / sinphi ) * sinphi
5343!
5344!-- Calculate potential visible diffuse radiation
5345    rdv = 0.4_wp * ( 600.0_wp - rdu ) * sinphi
5346!
5347!-- Calculate the water absorption in the-near infrared
5348    ww = 1320 * 10**( -1.195_wp + 0.4459_wp * log10( 1.0_wp / sinphi ) - 0.0345_wp * ( log10( 1.0_wp / sinphi ) )**2 )
5349!
5350!-- Calculate potential direct beam near-infrared radiation
5351    rdm = (720.0_wp * exp(-0.06_wp * (pres / pres_0) / sinphi ) - ww ) * sinphi     !< 720 = solar constant - 600
5352!
5353!-- Calculate potential diffuse near-infrared radiation
5354    rdn = 0.6_wp * ( 720 - rdm - ww ) * sinphi
5355!
5356!-- Compute visible and near-infrared radiation
5357    rv = MAX( 0.1_wp, rdu + rdv )
5358    rn = MAX( 0.01_wp, rdm + rdn )
5359!
5360!-- Compute ratio between input global radiation (here defined as solar radiation, dirict+diffuse)
5361!-- and total radiation computed here
5362    ratio = MIN( 0.89_wp, solar_rad / ( rv + rn ) )
5363!
5364!-- Calculate total visible radiation
5365    sv = ratio * rv
5366!
5367!-- Calculate fraction of par in the direct beam
5368    fv = MIN( 0.99_wp, ( 0.9_wp - ratio ) / 0.7_wp )              !< help variable
5369    fv = MAX( 0.01_wp, rdu / rv * ( 1.0_wp - fv**0.6667_wp ) )    !< fraction of par in the direct beam
5370!
5371!-- Compute direct and diffuse parts of par
5372    par_dir = fv * sv
5373    par_diff = sv - par_dir
5374
5375 END SUBROUTINE par_dir_diff
5376
5377 
5378 !-------------------------------------------------------------------
5379 !> rc_get_vpd: get vapour pressure deficit (kPa)
5380 !-------------------------------------------------------------------
5381 SUBROUTINE rc_get_vpd( temp, rh, vpd )
5382
5383!
5384!-- Input/output variables:
5385    REAL(wp), INTENT(IN) ::  temp    !< temperature (C)
5386    REAL(wp), INTENT(IN) ::  rh    !< relative humidity (%)
5387
5388    REAL(wp), INTENT(OUT) ::  vpd    !< vapour pressure deficit (kPa)
5389!
5390!-- Local variables:
5391    REAL(wp) ::  esat
5392!
5393!-- fit parameters:
5394    REAL(wp), PARAMETER ::  a1 = 6.113718e-01
5395    REAL(wp), PARAMETER ::  a2 = 4.43839e-02
5396    REAL(wp), PARAMETER ::  a3 = 1.39817e-03
5397    REAL(wp), PARAMETER ::  a4 = 2.9295e-05
5398    REAL(wp), PARAMETER ::  a5 = 2.16e-07
5399    REAL(wp), PARAMETER ::  a6 = 3.0e-09
5400!
5401!-- esat is saturation vapour pressure (kPa) at temp(C) following Monteith(1973)
5402    esat = a1 + a2 * temp + a3 * temp**2 + a4 * temp**3 + a5 * temp**4 + a6 * temp**5
5403    vpd  = esat * ( 1 - rh / 100 )
5404
5405 END SUBROUTINE rc_get_vpd
5406
5407
5408 !-------------------------------------------------------------------
5409 !> rc_gsoil_eff: compute effective soil conductance
5410 !-------------------------------------------------------------------
5411 SUBROUTINE rc_gsoil_eff( icmp, lu, sai, ust, nwet, t, gsoil_eff )
5412
5413!
5414!-- Input/output variables:
5415    INTEGER(iwp), INTENT(IN) ::  icmp          !< component index
5416    INTEGER(iwp), INTENT(IN) ::  lu            !< land use type, lu = 1,..., nlu
5417    INTEGER(iwp), INTENT(IN) ::  nwet          !< index for wetness
5418                                               !< nwet = 0 -> dry; nwet = 1 -> wet; nwet = 9 -> snow
5419                                               !< N.B. this routine cannot be called with nwet = 9,
5420                                               !< nwet = 9 should be handled outside this routine.
5421    REAL(wp), INTENT(IN) ::  sai               !< surface area index
5422    REAL(wp), INTENT(IN) ::  ust               !< friction velocity (m/s)
5423    REAL(wp), INTENT(IN) ::  t                 !< temperature (C)
5424    REAL(wp), INTENT(OUT) ::  gsoil_eff        !< effective soil conductance (m/s)
5425!
5426!-- local variables:
5427    REAL(wp) ::  rinc                          !< in canopy resistance  (s/m)
5428    REAL(wp) ::  rsoil_eff                     !< effective soil resistance (s/m)
5429!
5430!-- Soil resistance (numbers matched with lu_classes and component numbers)
5431    !     grs    ara    crp    cnf    dec    wat    urb   oth    des    ice    sav    trf    wai    med    sem
5432    REAL(wp), PARAMETER ::  rsoil(nlu_dep,ncmp) = reshape( (/ &
5433         1000.,  200.,  200.,  200.,  200., 2000.,  400., 1000., 2000., 2000., 1000.,  200., 2000.,  200.,  400., &    !< O3
5434         1000., 1000., 1000., 1000., 1000.,   10., 1000., 1000., 1000.,  500., 1000., 1000.,   10., 1000., 1000., &    !< SO2
5435         1000., 1000., 1000., 1000., 1000., 2000., 1000., 1000., 1000., 2000., 1000., 1000., 2000., 1000., 1000., &    !< NO2
5436         -999., -999., -999., -999., -999., 2000., 1000., -999., 2000., 2000., -999., -999., 2000., -999., -999., &    !< NO
5437         100.,  100.,  100.,  100.,  100.,   10.,  100.,  100.,  100., 1000.,  100.,  100.,   10.,  100.,  100.,  &    !< NH3
5438         -999., -999., -999., -999., -999., 2000., 1000.<