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

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

removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname

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