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

Last change on this file since 3664 was 3664, checked in by forkel, 3 years ago

Replaced misplaced location message by @todo

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