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

Last change on this file since 3796 was 3796, checked in by banzhafs, 3 years ago

Unused variables removed/taken care of from/in chemistry_model_mod

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