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

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

chemistry_model_mod code update to comply PALM coding rules

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