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

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

Removed forgotten write statements an some of the unused variables

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