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

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

read from unit 10 now also removed from subroutine chem_header

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