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

Last change on this file since 3655 was 3654, checked in by suehring, 3 years ago

Disable misplaced location message in chemistry_model_mod

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