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

Last change on this file since 3543 was 3543, checked in by suehring, 4 years ago

changes from last commit documented

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