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

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

clean-up debug prints

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