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

Last change on this file since 3767 was 3767, checked in by raasch, 3 years ago

unused variables removed from rrd-subroutines parameter list

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