Home

Context Navigation

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

Last change on this file since 3182 was 3182, checked in by suehring, 6 years ago
New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed
Property svn:keywords set to `Id`
File size: 76.3 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 1997-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	! Rename flags indicating outflow boundary conditions
25	!
26	! Former revisions:
27	! -----------------
28	! $Id: chemistry_model_mod.f90 3182 2018-07-27 13:36:03Z suehring $
29	! Revise output of surface quantities in case of overhanging structures
30	!
31	! 3045 2018-05-28 07:55:41Z Giersch
32	! error messages revised
33	!
34	! 3014 2018-05-09 08:42:38Z maronga
35	! Bugfix: nzb_do and nzt_do were not used for 3d data output
36	!
37	! 3004 2018-04-27 12:33:25Z Giersch
38	! Comment concerning averaged data output added
39	!
40	! 2932 2018-03-26 09:39:22Z maronga
41	! renamed chemistry_par to chemistry_parameters
42	!
43	! 2894 2018-03-15 09:17:58Z Giersch
44	! Calculations of the index range of the subdomain on file which overlaps with
45	! the current subdomain are already done in read_restart_data_mod,
46	! chem_last_actions was renamed to chem_wrd_local, chem_read_restart_data was
47	! renamed to chem_rrd_local, chem_write_var_list was renamed to
48	! chem_wrd_global, chem_read_var_list was renamed to chem_rrd_global,
49	! chem_skip_var_list has been removed, variable named found has been
50	! introduced for checking if restart data was found, reading of restart strings
51	! has been moved completely to read_restart_data_mod, chem_rrd_local is already
52	! inside the overlap loop programmed in read_restart_data_mod, todo list has
53	! bees extended, redundant characters in chem_wrd_local have been removed,
54	! the marker * end chemistry * is not necessary anymore, strings and their
55	! respective lengths are written out and read now in case of restart runs to
56	! get rid of prescribed character lengths
57	!
58	! 2815 2018-02-19 11:29:57Z suehring
59	! Bugfix in restart mechanism,
60	! rename chem_tendency to chem_prognostic_equations,
61	! implement vector-optimized version of chem_prognostic_equations,
62	! some clean up (incl. todo list)
63	!
64	! 2773 2018-01-30 14:12:54Z suehring
65	! Declare variables required for nesting as public
66	!
67	! 2772 2018-01-29 13:10:35Z suehring
68	! Bugfix in string handling
69	!
70	! 2768 2018-01-24 15:38:29Z kanani
71	! Shorten lines to maximum length of 132 characters
72	!
73	! 2766 2018-01-22 17:17:47Z kanani
74	! Removed preprocessor directive __chem
75	!
76	! 2756 2018-01-16 18:11:14Z suehring
77	! Fill values in 3D output introduced.
78	!
79	! 2718 2018-01-02 08:49:38Z maronga
80	! Initial revision
81	!
82	!
83	!
84	!
85	! Authors:
86	! --------
87	! @author Renate Forkel
88	! @author Farah Kanani-Suehring
89	! @author Klaus Ketelsen
90	! @author Basit Khan
91	!
92	!
93	!------------------------------------------------------------------------------!
94	! Description:
95	! ------------
96	!> Chemistry model for PALM-4U
97	!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
98	!> allowed to use the chemistry model in a precursor run and additionally
99	!> not using it in a main run
100	!> @todo Update/clean-up todo list! (FK)
101	!> @todo Set proper fill values (/= 0) for chem output arrays! (FK)
102	!> @todo Add routine chem_check_parameters, add checks for inconsistent or
103	!> unallowed parameter settings.
104	!> CALL of chem_check_parameters from check_parameters. (FK)
105	!> @todo Make routine chem_header available, CALL from header.f90
106	!> (see e.g. how it is done in routine lsm_header in
107	!> land_surface_model_mod.f90). chem_header should include all setup
108	!> info about chemistry parameter settings. (FK)
109	!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
110	!> @todo Separate boundary conditions for each chem spcs to be implemented
111	!> @todo Currently only total concentration are calculated. Resolved, parameterized
112	!> and chemistry fluxes although partially and some completely coded but
113	!> are not operational/activated in this version. bK.
114	!> @todo slight differences in passive scalar and chem spcs when chem reactions
115	!> turned off. Need to be fixed. bK
116	!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
117	!> @todo subroutine set_const_initial_values to be taken out from chemistry_model_mod !bK.
118	!> @todo chemistry error messages
119	!> @todo Format this module according to PALM coding standard (see e.g. module
120	!> template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or
121	!> D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16)
122	!
123	!------------------------------------------------------------------------------!
124
125	MODULE chemistry_model_mod
126
127	USE kinds, ONLY: wp, iwp
128	USE indices, ONLY: nz, nzb,nzt,nysg,nyng,nxlg,nxrg,nys,nyn
129	USE pegrid, ONLY: myid, threads_per_task
130	USE control_parameters, ONLY: dt_3d, ws_scheme_sca, initializing_actions, message_string, &
131	omega, tsc, intermediate_timestep_count, &
132	intermediate_timestep_count_max, &
133	timestep_scheme, use_prescribed_profile_data
134	USE arrays_3d, ONLY: hyp, pt, rdf_sc, tend, zu
135	USE chem_gasphase_mod, ONLY: nspec, spc_names, nkppctrl, nmaxfixsteps, &
136	t_steps, fill_temp, chem_gasphase_integrate, &
137	nvar, atol, rtol, nphot, phot_names
138	USE cpulog, ONLY: cpu_log, log_point
139
140	USE chem_modules
141
142
143	IMPLICIT NONE
144	PRIVATE
145	SAVE
146
147	!- Define chemical variables
148
149	TYPE species_def
150	CHARACTER(LEN=8) :: name
151	CHARACTER(LEN=16) :: unit
152	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc
153	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_av
154	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_p
155	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: tconc_m
156	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: cssws_av
157	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: flux_s_cs
158	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: diss_s_cs
159	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: flux_l_cs
160	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: diss_l_cs
161	REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: conc_pr_init
162	END TYPE species_def
163
164
165	TYPE photols_def
166	CHARACTER(LEN=8) :: name
167	CHARACTER(LEN=16) :: unit
168	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: freq
169	END TYPE photols_def
170
171
172	PUBLIC species_def
173	PUBLIC photols_def
174
175
176	TYPE(species_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: chem_species
177	TYPE(photols_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: phot_frequen
178
179	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_1
180	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_2
181	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_3
182	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_av
183	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: freq_1
184
185	INTEGER,DIMENSION(nkppctrl) :: icntrl ! Fine tuning kpp
186	REAL(kind=wp),DIMENSION(nkppctrl) :: rcntrl ! Fine tuning kpp
187
188	CHARACTER(10), PUBLIC :: photolysis_scheme
189	! 'constant',
190	! 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180
191	! 'fastj' (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
192
193	PUBLIC nspec
194	PUBLIC nvar
195	PUBLIC spc_names
196	PUBLIC spec_conc_2
197
198	!- Interface section
199	INTERFACE chem_boundary_conds
200	MODULE PROCEDURE chem_boundary_conds
201	END INTERFACE chem_boundary_conds
202
203	INTERFACE chem_init
204	MODULE PROCEDURE chem_init
205	END INTERFACE chem_init
206
207	INTERFACE chem_init_profiles
208	MODULE PROCEDURE chem_init_profiles
209	END INTERFACE chem_init_profiles
210
211	INTERFACE chem_parin
212	MODULE PROCEDURE chem_parin
213	END INTERFACE chem_parin
214
215	INTERFACE chem_integrate
216	MODULE PROCEDURE chem_integrate_ij
217	END INTERFACE chem_integrate
218
219	INTERFACE chem_swap_timelevel
220	MODULE PROCEDURE chem_swap_timelevel
221	END INTERFACE chem_swap_timelevel
222
223	INTERFACE chem_define_netcdf_grid
224	MODULE PROCEDURE chem_define_netcdf_grid
225	END INTERFACE chem_define_netcdf_grid
226
227	INTERFACE chem_data_output_3d
228	MODULE PROCEDURE chem_data_output_3d
229	END INTERFACE chem_data_output_3d
230
231	INTERFACE chem_check_data_output
232	MODULE PROCEDURE chem_check_data_output
233	END INTERFACE chem_check_data_output
234
235	INTERFACE chem_check_data_output_pr
236	MODULE PROCEDURE chem_check_data_output_pr
237	END INTERFACE chem_check_data_output_pr
238
239	INTERFACE chem_3d_data_averaging
240	MODULE PROCEDURE chem_3d_data_averaging
241	END INTERFACE chem_3d_data_averaging
242
243	INTERFACE chem_wrd_local
244	MODULE PROCEDURE chem_wrd_local
245	END INTERFACE chem_wrd_local
246
247	INTERFACE chem_rrd_local
248	MODULE PROCEDURE chem_rrd_local
249	END INTERFACE chem_rrd_local
250
251	INTERFACE chem_prognostic_equations
252	MODULE PROCEDURE chem_prognostic_equations
253	MODULE PROCEDURE chem_prognostic_equations_ij
254	END INTERFACE chem_prognostic_equations
255
256	INTERFACE chem_header
257	MODULE PROCEDURE chem_header
258	END INTERFACE chem_header
259
260	INTERFACE chem_emissions
261	MODULE PROCEDURE chem_emissions
262	END INTERFACE chem_emissions
263
264	! INTERFACE chem_wrd_global
265	! MODULE PROCEDURE chem_wrd_global
266	! END INTERFACE chem_wrd_global
267	!
268	! INTERFACE chem_rrd_global
269	! MODULE PROCEDURE chem_rrd_global
270	! END INTERFACE chem_rrd_global
271
272
273	PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
274	chem_check_data_output_pr, chem_data_output_3d, &
275	chem_define_netcdf_grid, chem_emissions, chem_header, chem_init, &
276	chem_init_profiles, chem_integrate, chem_wrd_local, &
277	chem_parin, chem_prognostic_equations, &
278	chem_rrd_local, chem_swap_timelevel
279
280
281	CONTAINS
282
283	!------------------------------------------------------------------------------!
284	!
285	! Description:
286	! ------------
287	!> Subroutine to initialize and set all boundary conditions for chemical species
288	!------------------------------------------------------------------------------!
289
290	SUBROUTINE chem_boundary_conds( mode )
291
292	USE control_parameters, &
293	ONLY: air_chemistry, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
294	bc_radiation_s
295	USE indices, &
296	ONLY: nxl, nxr, nxlg, nxrg, nyng, nysg, nzt
297
298	! USE prognostic_equations_mod, &
299
300	USE arrays_3d, &
301	ONLY: dzu
302	USE surface_mod, &
303	ONLY: bc_h
304
305	CHARACTER (len=*), INTENT(IN) :: mode
306	INTEGER(iwp) :: i !< grid index x direction.
307	INTEGER(iwp) :: j !< grid index y direction.
308	INTEGER(iwp) :: k !< grid index z direction.
309	INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing.
310	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls.
311	INTEGER(iwp) :: m !< running index surface elements.
312	INTEGER(iwp) :: lsp !< running index for chem spcs.
313	INTEGER(iwp) :: lph !< running index for photolysis frequencies
314
315
316	SELECT CASE ( TRIM( mode ) )
317	CASE ( 'init' )
318	DO lsp = 1, nspec
319	IF ( surface_csflux(lsp) == 9999999.9_wp ) THEN
320	constant_csflux(lsp) = .FALSE.
321	ENDIF
322	ENDDO
323
324	IF ( bc_cs_b == 'dirichlet' ) THEN
325	ibc_cs_b = 0
326	ELSEIF ( bc_cs_b == 'neumann' ) THEN
327	ibc_cs_b = 1
328	ELSE
329	! message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' ! bK commented
330	CALL message( 'chem_boundary_conds', 'CM0010', 1, 2, 0, 6, 0 ) !< chemistry_model_mod should have special error numbers --> "CHEM###",
331	ENDIF
332	!
333	!-- Set Integer flags and check for possible erroneous settings for top
334	!-- boundary condition. bK added _cs_ here.
335	IF ( bc_cs_t == 'dirichlet' ) THEN
336	ibc_cs_t = 0
337	ELSEIF ( bc_cs_t == 'neumann' ) THEN
338	ibc_cs_t = 1
339	ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN
340	ibc_cs_t = 2
341	ELSEIF ( bc_cs_t == 'nested' ) THEN
342	ibc_cs_t = 3
343	ELSE
344	! message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'
345	CALL message( 'check_parameters', 'CM0011', 1, 2, 0, 6, 0 )
346	ENDIF
347
348
349	CASE ( 'set_bc_bottomtop' )
350	!-- Bottom boundary condtions for chemical species
351	DO lsp = 1, nspec
352	IF ( ibc_cs_b == 0 ) THEN
353	DO l = 0, 1
354	!-- Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
355	!-- the chem_species(nsp)%conc_p value at the topography top (k-1)
356	!-- is set; for downward-facing surfaces (l=1), kb=1, i.e. the
357	!-- value at the topography bottom (k+1) is set.
358
359	kb = MERGE( -1, 1, l == 0 )
360	!$OMP PARALLEL DO PRIVATE( i, j, k )
361	DO m = 1, bc_h(l)%ns
362	i = bc_h(l)%i(m)
363	j = bc_h(l)%j(m)
364	k = bc_h(l)%k(m)
365	chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i)
366	ENDDO
367	ENDDO
368
369	ELSEIF ( ibc_cs_b == 1 ) THEN
370	! in boundary_conds there is som extra loop over m here for passive tracer
371	DO l = 0, 1
372	kb = MERGE( -1, 1, l == 0 )
373	!$OMP PARALLEL DO PRIVATE( i, j, k )
374	DO m = 1, bc_h(l)%ns
375	i = bc_h(l)%i(m)
376	j = bc_h(l)%j(m)
377	k = bc_h(l)%k(m)
378	chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i)
379
380	!< @todo: chem_species(nsp)%conc_p(k+kb,j,i) = chem_species(nsp)%conc(k+kb,j,i), &
381	!< pls loop over nsp=1, NSPEC.
382	!< @todo: We should also think about the possibility to have &
383	!< individual boundary conditions for each species? This means, bc_cs_b, &
384	!< bc_cs_t, ibc_cs_b, ibc_cs_t would need to be added to TYPE chem_species(nsp)%, &
385	!< and the loop over nsp would have to be put outside of this IF-clause.
386	!< i think its better we keep the same bonundary cond i.e. dirichlet or neumann
387	!< for all chem spcs. ... !bK
388
389	ENDDO
390	ENDDO
391	ENDIF
392	ENDDO ! end lsp loop
393
394	!-- Top boundary conditions for chemical species - Should this not be done for all species?
395	IF ( ibc_cs_t == 0 ) THEN
396	DO lsp = 1, nspec
397	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)
398	ENDDO
399	ELSEIF ( ibc_cs_t == 1 ) THEN
400	DO lsp = 1, nspec
401	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
402	ENDDO
403	ELSEIF ( ibc_cs_t == 2 ) THEN
404	DO lsp = 1, nspec
405	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
406	ENDDO
407	!@todo: bc_cs_t_val needs to be calculated,
408	!see bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
409	!(in time_integration). pt_init would be the counterpart to
410	!chem_species(i)%conc_pr_init (see kchem_driver_FKa1408.f90 of my
411	!"Hints: prescribing initial concentration.
412	ENDIF
413	!
414	CASE ( 'set_bc_lateral' ) ! bK commented it
415	!-- Lateral boundary conditions for chem species at inflow boundary
416	!-- are automatically set when chem_species concentration is
417	!-- initialized. The initially set value at the inflow boundary is not
418	!-- touched during time integration, hence, this boundary value remains
419	!-- at a constant value, which is the concentration that flows into the
420	!-- domain.
421	!-- Lateral boundary conditions for chem species at outflow boundary
422
423	IF ( bc_radiation_s ) THEN
424	DO lsp = 1, nspec
425	chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
426	ENDDO
427	ELSEIF ( bc_radiation_n ) THEN
428	DO lsp = 1, nspec
429	chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
430	ENDDO
431	ELSEIF ( bc_radiation_l ) THEN
432	DO lsp = 1, nspec
433	chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
434	ENDDO
435	ELSEIF ( bc_radiation_r ) THEN
436	DO lsp = 1, nspec
437	chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
438	ENDDO
439	ENDIF
440
441	END SELECT
442
443	END SUBROUTINE chem_boundary_conds
444	!
445	!------------------------------------------------------------------------------!
446	!
447	! Description:
448	! ------------
449	!> Subroutine defining initial vertical profiles of chemical species (given by
450	!> namelist parameters chem_profiles and chem_heights) --> which should work
451	!> analogue to parameters u_profile, v_profile and uv_heights)
452	!------------------------------------------------------------------------------!
453	SUBROUTINE chem_init_profiles !< SUBROUTINE is called from chem_init in case of
454	!< TRIM( initializing_actions ) /= 'read_restart_data'
455	!< We still need to see what has to be done in case of restart run
456	USE chem_modules
457
458	IMPLICIT NONE
459
460	!-- Local variables
461	INTEGER :: lsp !< running index for number of species in derived data type species_def
462	INTEGER :: lsp_pr !< running index for number of species in cs_names, cs_profiles etc
463	INTEGER :: lpr_lev !< running index for profile level for each chem spcs.
464	INTEGER :: npr_lev !< the next available profile lev
465
466	!-----------------
467	!-- To prescribe/initialize a vertically constant 'cs_profile', another parameter
468	!-- "cs_surface" is introduced. If "cs_profile" and "cs_heights" are prescribed,
469	!-- their values will override the constant profile given by "cs_surface".
470
471	! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
472	lsp_pr = 1
473	DO WHILE ( TRIM( cs_name( lsp_pr ) ) /= 'novalue' ) !'novalue' is the default
474	DO lsp = 1, nspec !
475	!-- Create initial profile (conc_pr_init) for each chemical species
476	IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_pr) ) ) THEN !
477	! IF ( cs_profile(1,1) == 9999999.9_wp ) THEN
478	!-- Set a vertically constant profile based on the surface conc (cs_surface(lsp_pr)) of each species
479	DO lpr_lev = 0, nzt+1
480	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_pr)
481	ENDDO
482
483	! ELSE
484	! IF ( cs_heights(lsp,1) /= 0.0_wp ) THEN
485	! message_string = 'cs_heights(1,1) must be 0.0'
486	! CALL message( 'chem_check_parameters', 'CM0012', 1, 2, 0, 6, 0 )
487	! ENDIF
488	!
489	! IF ( omega /= 0.0_wp ) THEN
490	! message_string = 'Coriolis force must be switched off (by setting omega=0.0)' // &
491	! ' when prescribing the forcing by u_profile and v_profile'
492	! CALL message( 'check_parameters', 'PA0347', 1, 2, 0, 6, 0 )
493	! ENDIF
494	!
495	! use_prescribed_profile_data = .TRUE.
496	!
497	! npr_lev = 1
498	! ! chem_sddpecies(lsp)%conc_pr_init(0) = 0.0_wp
499	! DO lpr_lev = 1, nz+1
500	! IF ( npr_lev < 100 ) THEN
501	! DO WHILE ( cs_heights(lsp, npr_lev+1) <= zu(lpr_lev) )
502	! npr_lev = npr_lev + 1
503	! IF ( npr_lev == 100 ) EXIT
504	! ENDDO
505	! ENDIF
506	!
507	! IF ( npr_lev < 100 .AND. cs_heights(lsp, npr_lev + 1) /= 9999999.9_wp ) THEN
508	! chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp, npr_lev) + &
509	! ( zu(lpr_lev) - cs_heights(lsp, npr_lev) ) / &
510	! ( cs_heights(lsp, npr_lev + 1) - cs_heights(lsp, npr_lev ) ) * &
511	! ( cs_profile(lsp, npr_lev + 1) - cs_profile(lsp, npr_lev ) )
512	! ELSE
513	! chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp, npr_lev)
514	! ENDIF
515	! ENDDO
516	! ENDIF
517	!-- If a profile is prescribed explicity using cs_profiles and cs_heights, we then have to fill the
518	!-- chem_species(lsp)%conc_pr_init for the specific "lsp" based on the cs_profiles(lsp_pr,:)
519	!-- and cs_heights(lsp_pr,:).
520	ENDIF
521	ENDDO
522	lsp_pr = lsp_pr + 1
523	ENDDO
524	! ELSE
525	! IF (chem_debug1 ) print*,'code to be added for initializing_actions == read_restart_data' !bK
526	! ENDIF
527
528	!-- Now, go back to chem_init and use the contents of chem_species(lsp)%conc_pr_init to
529	!-- initialize the 3D conc arrays, as it is currently taken care of in chem_set_constant_values.
530	!-- After initializing the 3D arrays, these can be used to set the boundary conditions in the
531	!-- subroutine kchem_boundary_conds, which should be called from boundary_conds.f90.
532
533
534	END SUBROUTINE chem_init_profiles
535	!
536	!------------------------------------------------------------------------------!
537	!
538	! Description:
539	! ------------
540	!> Subroutine initializating chemistry_model_mod
541	!------------------------------------------------------------------------------!
542	SUBROUTINE chem_init
543
544
545	USE control_parameters, &
546	ONLY: message_string, io_blocks, io_group, turbulent_inflow
547
548	USE arrays_3d, &
549	ONLY: mean_inflow_profiles
550
551	USE pegrid
552
553	IMPLICIT none
554	!-- local variables
555	INTEGER :: i,j !< running index for for horiz numerical grid points
556	INTEGER :: lsp !< running index for chem spcs
557	INTEGER :: lpr_lev !< running index for chem spcs profile level
558	!
559	!-- NOPOINTER version not implemented yet
560	! #if defined( __nopointer )
561	! message_string = 'The chemistry module only runs with POINTER version'
562	! CALL message( 'chemistry_model_mod', 'CM0001', 1, 2, 0, 6, 0 )
563	! #endif
564	!
565	!-- Allocate memory for chemical species
566	ALLOCATE( chem_species(nspec) )
567	ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
568	ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
569	ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
570	ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
571	ALLOCATE( phot_frequen(nphot) )
572	ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
573	ALLOCATE( bc_cs_t_val(nspec) )
574	!
575	!-- Initialize arrays
576	spec_conc_1 (:,:,:,:) = 0.0_wp
577	spec_conc_2 (:,:,:,:) = 0.0_wp
578	spec_conc_3 (:,:,:,:) = 0.0_wp
579	spec_conc_av(:,:,:,:) = 0.0_wp
580
581
582	DO lsp = 1, nspec
583	chem_species(lsp)%name = spc_names(lsp)
584
585	chem_species(lsp)%conc (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1 (:,:,:,lsp)
586	chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2 (:,:,:,lsp)
587	chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_3 (:,:,:,lsp)
588	chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_av(:,:,:,lsp)
589
590	ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))
591	chem_species(lsp)%cssws_av = 0.0_wp
592	!
593	!-- (todo (FK): This needs to be revised. This block must go somewhere else)
594	! IF ( ws_scheme_sca ) THEN
595	ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))
596	ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))
597	ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
598	ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
599	chem_species(lsp)%flux_s_cs = 0.0_wp
600	chem_species(lsp)%flux_l_cs = 0.0_wp
601	chem_species(lsp)%diss_s_cs = 0.0_wp
602	chem_species(lsp)%diss_l_cs = 0.0_wp
603	! ENDIF
604	!
605	!-- Allocate memory for initial concentration profiles
606	!-- (concentration values come from namelist)
607	!-- (todo (FK): Because of this, chem_init is called in palm before
608	!-- check_parameters, since conc_pr_init is used there.
609	!-- We have to find another solution since chem_init should
610	!-- eventually be called from init_3d_model!!)
611	ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
612	chem_species(lsp)%conc_pr_init(:) = 0.0_wp
613
614	ENDDO
615
616	!
617	!-- Set initial concentration of profiles prescribed by parameters cs_profile
618	!-- and cs_heights in the namelist &chemistry_parameters
619	!-- (todo (FK): chem_init_profiles not ready yet, has some bugs)
620	! CALL chem_init_profiles
621	!
622	!-- Initialize model variables
623
624
625	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
626	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
627
628
629	!-- First model run of a possible job queue.
630	!-- Initial profiles of the variables must be computes.
631	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
632	! CALL location_message( 'initializing with 1D model profiles', .FALSE. )
633	!
634	! CALL init_1d_model ...... decide to call or not later !bK
635
636	!-- Transfer initial profiles to the arrays of the 3D model
637	DO lsp = 1, nspec
638	DO i = nxlg, nxrg
639	DO j = nysg, nyng
640	DO lpr_lev = 1, nz + 1
641	chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
642	ENDDO
643	ENDDO
644	ENDDO
645	ENDDO
646
647	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
648	THEN
649	! CALL location_message( 'initializing with constant profiles', .FALSE. )
650
651
652
653	!-- Set initial horizontal velocities at the lowest computational grid
654	!-- levels to zero in order to avoid too small time steps caused by the
655	!-- diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
656	!-- in the limiting formula!).
657
658	DO lsp = 1, nspec
659	DO i = nxlg, nxrg
660	DO j = nysg, nyng
661	chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init !ITS THERE bK
662	ENDDO
663	ENDDO
664	ENDDO
665
666	! ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
667	! THEN
668	! CALL location_message( 'initializing by user', .FALSE. )
669	!
670	!-- Initialization will completely be done by the user
671	!-- (FK: This should be called only once, in init_3d_model, i.e. remove it here)
672	! CALL user_init_3d_model
673	! CALL location_message( 'finished', .TRUE. )
674
675	ENDIF
676
677	!-- Store initial chem spcs profile
678	! DO lsp = 1, nvar
679	! hom_cs(:,1,115,:) = SPREAD( chem_species(lsp)%conc(:,nys,nxl), 2, statistic_regions+1 )
680	! ENDDO
681	!
682	!-- If required, change the surface chem spcs at the start of the 3D run
683	IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN
684	DO lsp = 1, nspec
685	chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) + &
686	cs_surface_initial_change(lsp)
687	ENDDO
688	ENDIF
689	!
690	!-- Initiale old and new time levels.
691	DO lsp = 1, nvar
692	chem_species(lsp)%tconc_m = 0.0_wp
693	chem_species(lsp)%conc_p = chem_species(lsp)%conc
694	ENDDO
695
696	ENDIF
697
698
699
700	!--- new code add above this line
701	DO lsp = 1, nphot
702	phot_frequen(lsp)%name = phot_names(lsp)
703	! IF( myid == 0 ) THEN
704	! WRITE(6,'(a,i4,3x,a)') 'Photolysis: ',lsp,trim(phot_names(lsp))
705	! ENDIF
706	phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => freq_1(:,:,:,lsp)
707	ENDDO
708
709	!-- Set initial values
710	! Not required any more, this can now be done with the namelist by setting cs_surface
711	! and cs_name without specifying cs_profile (Nevertheless, we still want to keep it for a while)
712	! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
713	! CALL set_const_initial_values
714	! ENDIF
715
716	RETURN
717
718	CONTAINS
719	!------------------------------------------------------------------------------!
720	!
721	! Description:
722	! ------------
723	!> Subroutine setting constant initial values of chemical species
724	!------------------------------------------------------------------------------!
725	SUBROUTINE set_const_initial_values
726	! Not required any more, this can now be done with the namelist by setting cs_surface
727	! and cs_name without specifying cs_profile (Nevertheless, we still want to keep it for a while)
728	IMPLICIT none
729
730	!-- local variables
731	INTEGER :: lsp
732
733	IF(myid == 0) THEN
734	write(6,'(/,a,/)') ' chemics >>>> Set constant Initial Values: '
735	ENDIF
736
737	! Default values are taken from smog.def from supplied kpp example
738	DO lsp = 1, nspec
739	IF(trim(chem_species(lsp)%name) == 'NO') THEN
740	! chem_species(lsp)%conc = 8.725*1.0E+08
741	! chem_species(lsp)%conc = 0.05_wp !added by bK
742	chem_species(lsp)%conc = 0.01_wp !added by RFo
743	ELSEIF (trim(chem_species(lsp)%name) == 'NO2') THEN
744	! chem_species(lsp)%conc = 2.240*1.0E+08
745	! chem_species(lsp)%conc = 0.01_wp !added by bK
746	chem_species(lsp)%conc = 0.05_wp !added by RFo
747	ELSEIF( trim( chem_species(lsp)%name ) == 'O3' ) THEN
748	chem_species(lsp)%conc = 0.05_wp !added by bK
749	ELSEIF(trim(chem_species(lsp)%name) == 'H2O') THEN
750	! chem_species(lsp)%conc = 5.326*1.0E+11
751	chem_species(lsp)%conc = 1.30*1.0E+4_wp !added by bK
752	ELSEIF(trim(chem_species(lsp)%name) == 'O2') THEN
753	chem_species(lsp)%conc = 2.0*1.0E+5_wp !added by bK
754	ELSEIF(trim(chem_species(lsp)%name) == 'RH') THEN
755	chem_species(lsp)%conc = 0.001_wp !added by RFo
756	ELSEIF(trim(chem_species(lsp)%name) == 'CO') THEN
757	chem_species(lsp)%conc = 0.5_wp !added by RFo
758	ELSEIF(trim(chem_species(lsp)%name) == 'RCHO') THEN
759	! chem_species(lsp)%conc = 2.0_wp !added by bK
760	chem_species(lsp)%conc = 0.01_wp !added by RFo
761	! ELSEIF(trim(chem_species(lsp)%name) == 'OH') THEN
762	! chem_species(lsp)%conc = 1.0*1.0E-07_wp !added by bK
763	! ELSEIF(trim(chem_species(lsp)%name) == 'HO2') THEN
764	! chem_species(lsp)%conc = 1*1.0E-7_wp !added by bK
765	! ELSEIF(trim(chem_species(lsp)%name) == 'RCOO2') THEN ! corrected RFo
766	! chem_species(lsp)%conc = 1.0*1.0E-7_wp !added by bK
767	! ELSEIF(trim(chem_species(lsp)%name) == 'RCOO2NO2') THEN
768	! chem_species(lsp)%conc = 1.0*1.0E-7_wp !added by bK
769	ELSE
770	! H2O = 2.0e+04;
771	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) = 0.0_wp
772	ENDIF
773
774	IF(myid == 0) THEN
775	WRITE(6,'(a,3x,a,3x,a,e12.4)') ' Species: ',chem_species(lsp)%name(1:7), &
776	'Initial Value = ',chem_species(lsp)%conc(nzb,nysg,nxlg)
777	ENDIF
778	ENDDO
779
780	! #if defined( __nopointer )
781	! !kk Hier mit message abbrechen
782	! if(myid == 0) then
783	! write(6,*) ' KPP does only run with POINTER Version'
784	! end if
785	! stop 'error'
786	! #endif
787
788	RETURN
789	END SUBROUTINE set_const_initial_values
790
791
792	END SUBROUTINE chem_init
793	!
794	!------------------------------------------------------------------------------!
795	!
796	! Description:
797	! ------------
798	!> Subroutine defining parin for &chemistry_parameters for chemistry model
799	!------------------------------------------------------------------------------!
800	SUBROUTINE chem_parin
801
802	USE control_parameters, &
803	ONLY: air_chemistry
804
805	USE chem_modules
806
807	USE kinds
808
809	IMPLICIT none
810
811	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
812
813	REAL(wp), DIMENSION(nmaxfixsteps) :: my_steps !< List of fixed timesteps my_step(1) = 0.0 automatic stepping
814
815	NAMELIST /chemistry_parameters/ bc_cs_b, &
816	bc_cs_t, &
817	call_chem_at_all_substeps, &
818	chem_debug0, &
819	chem_debug1, &
820	chem_debug2, &
821	chem_gasphase_on, &
822	cs_heights, &
823	cs_name, &
824	cs_profile, &
825	cs_profile_name, &
826	cs_surface, &
827	emiss_factor_main, &
828	emiss_factor_side, &
829	icntrl, &
830	main_street_id, &
831	max_street_id, &
832	my_steps, &
833	rcntrl, &
834	side_street_id, &
835	photolysis_scheme, &
836	wall_csflux, &
837	cs_vertical_gradient, &
838	top_csflux, &
839	surface_csflux, &
840	surface_csflux_name, &
841	cs_surface_initial_change, &
842	cs_vertical_gradient_level
843
844	!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
845	!-- so this way we could prescribe a specific flux value for each species
846	!> chemistry_parameters for initial profiles
847	!> cs_names = 'O3', 'NO2', 'NO', ... to set initial profiles)
848	!> cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
849	!> cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
850	!> cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, ..... (chem spcs conc at height lvls chem_heights(1,:)) etc.
851	!> If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
852	!> then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
853
854	!
855	!-- Read chem namelist
856	!-- (todo: initialize these parameters in declaration part, do this for
857	!-- all chemistry_parameters namelist parameters)
858	icntrl = 0
859	rcntrl = 0.0_wp
860	my_steps = 0.0_wp
861	icntrl(2) = 1
862	photolysis_scheme = 'simple'
863	atol = 1.0_wp
864	rtol = 0.01_wp
865	!
866	!-- Try to find chemistry package
867	REWIND ( 11 )
868	line = ' '
869	DO WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
870	READ ( 11, '(A)', END=10 ) line
871	ENDDO
872	BACKSPACE ( 11 )
873	!
874	!-- Read chemistry namelist
875	READ ( 11, chemistry_parameters )
876	!
877	!-- Enable chemistry model
878	air_chemistry = .TRUE.
879
880
881	10 CONTINUE
882
883	t_steps = my_steps !(todo: Why not directly make t_steps a
884	! namelist parameter?)
885
886	END SUBROUTINE chem_parin
887
888	!
889	!------------------------------------------------------------------------------!
890	!
891	! Description:
892	! ------------
893	!> Subroutine to integrate chemical species in the given chemical mechanism
894	!------------------------------------------------------------------------------!
895
896	SUBROUTINE chem_integrate_ij (i, j)
897
898	USE cpulog, &
899	ONLY: cpu_log, log_point
900	USE statistics, & ! ## RFo
901	ONLY: weight_pres
902	USE control_parameters, & ! ## RFo
903	ONLY: dt_3d, intermediate_timestep_count
904
905	IMPLICIT none
906	INTEGER,INTENT(IN) :: i,j
907
908	!-- local variables
909	INTEGER(iwp) :: lsp !< running index for chem spcs.
910	INTEGER(iwp) :: lph !< running index for photolysis frequencies
911	INTEGER :: k,m,istatf
912	INTEGER,dimension(20) :: istatus
913	REAL(kind=wp),dimension(nzb+1:nzt,nspec) :: tmp_conc
914	REAL(kind=wp),dimension(nzb+1:nzt) :: tmp_temp
915	REAL(kind=wp),dimension(nzb+1:nzt,nphot) :: tmp_phot
916	REAL(kind=wp),dimension(nzb+1:nzt) :: tmp_fact
917	REAL(kind=wp),dimension(nzb+1:nzt) :: tmp_fact_i !< conversion factor between molecules cm^{-3} and ppm
918	REAL(wp) :: conv !< conversion factor
919	REAL(wp), PARAMETER :: ppm2fr = 1.0e-6_wp !< Conversion factor ppm to fraction
920	REAL(wp), PARAMETER :: pref_i = 1._wp / 100000.0_wp !< inverse reference pressure (1/Pa)
921	REAL(wp), PARAMETER :: t_std = 273.15_wp !< standard pressure (Pa)
922	REAL(wp), PARAMETER :: p_std = 101325.0_wp !< standard pressure (Pa)
923	REAL(wp), PARAMETER :: r_cp = 0.286_wp !< R / cp (exponent for potential temperature)
924	REAL(wp), PARAMETER :: vmolcm = 22.414e3_wp !< Mole volume (22.414 l) in cm^{-3}
925	REAL(wp), PARAMETER :: xna = 6.022e23_wp !< Avogadro number (molecules/mol)
926
927
928	REAL(kind=wp) :: dt_chem
929
930	CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' )
931	!< Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
932	IF (chem_gasphase_on) THEN
933
934	tmp_temp(:) = pt(:,j,i) * ( hyp(nzb+1:nzt) / 100000.0_wp )**0.286_wp
935	! ppm to molecules/cm**3
936	! tmp_fact = 1.e-6_wp6.022e23_wp/(22.414_wp1000._wp) * 273.15_wp * hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp )
937	conv = ppm2fr * xna / vmolcm
938	tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
939	tmp_fact_i = 1.0_wp/tmp_fact
940
941	CALL fill_temp (istatf, tmp_temp) !< Load constant temperature into kpp context
942	! CALL fill_temp (istatf, pt(nzb+1:nzt,j,i)) !< Load temperature into kpp context
943
944	DO lsp = 1,nspec
945	tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
946	ENDDO
947
948	DO lph = 1,nphot
949	tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)
950	ENDDO
951
952	IF(myid == 0 .AND. chem_debug0 ) THEN
953	IF (i == 10 .and. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
954	ENDIF
955
956	!-- Compute length of time step # RFo
957	IF ( call_chem_at_all_substeps ) THEN
958	dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
959	ELSE
960	dt_chem = dt_3d
961	ENDIF
962
963	CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' )
964
965
966	CALL chem_gasphase_integrate (dt_chem, tmp_conc, tmp_temp, tmp_phot, istatus=istatus)
967
968
969	CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' )
970
971	DO lsp = 1,nspec
972	chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:) ! RFo
973	ENDDO
974
975	! IF (myid == 0 .AND. chem_debug2 ) THEN
976	! IF (i == 10 .and. j == 10) WRITE(6,'(a,8i7)') ' KPP Status ',istatus(1:8)
977	! write(6,'(a,8i7)') ' KPP Status ',istatus(1:8)
978	! ENDIF
979
980	ENDIF
981	CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' )
982
983	RETURN
984	END SUBROUTINE chem_integrate_ij
985	!
986	!------------------------------------------------------------------------------!
987	!
988	! Description:
989	! ------------
990	!> Subroutine for swapping of timelevels for chemical species
991	!> called out from subroutine swap_timelevel
992	!------------------------------------------------------------------------------!
993
994	SUBROUTINE chem_swap_timelevel (level)
995	IMPLICIT none
996
997	INTEGER,INTENT(IN) :: level
998
999	!-- local variables
1000
1001	INTEGER :: lsp
1002
1003	! print,' entering chem_swap_timelevel *) '
1004	if(level == 0) then
1005	do lsp=1, nvar
1006	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp)
1007	chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp)
1008	end do
1009	else
1010	do lsp=1, nvar
1011	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp)
1012	chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp)
1013	end do
1014	end if
1015
1016	RETURN
1017	END SUBROUTINE chem_swap_timelevel
1018
1019	!------------------------------------------------------------------------------!
1020	!
1021	! Description:
1022	! ------------
1023	!> Subroutine defining appropriate grid for netcdf variables.
1024	!> It is called out from subroutine netcdf.
1025	!------------------------------------------------------------------------------!
1026	SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1027
1028	IMPLICIT NONE
1029
1030	CHARACTER (LEN=*), INTENT(IN) :: var !<
1031	LOGICAL, INTENT(OUT) :: found !<
1032	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
1033	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
1034	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
1035
1036	found = .TRUE.
1037
1038	if(var(1:3) == 'kc_') then !< always the same grid for chemistry variables
1039	grid_x = 'x'
1040	grid_y = 'y'
1041	grid_z = 'zu' !< kk Use same z axis as u variables. Has to be checked if OK
1042	else
1043	found = .FALSE.
1044	grid_x = 'none'
1045	grid_y = 'none'
1046	grid_z = 'none'
1047	end if
1048
1049	! write(6,*) 'chem_define_netcdf_grid ',TRIM(var),' ',trim(grid_x),' ',found
1050
1051	END SUBROUTINE chem_define_netcdf_grid
1052	!
1053	!------------------------------------------------------------------------------!
1054	!
1055	! Description:
1056	! ------------
1057	!> Subroutine for checking data output for chemical species
1058	!------------------------------------------------------------------------------!
1059
1060	SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
1061
1062
1063	USE control_parameters, &
1064	ONLY: data_output, message_string
1065
1066	IMPLICIT NONE
1067
1068	CHARACTER (LEN=*) :: unit !<
1069	CHARACTER (LEN=*) :: var !<
1070
1071	INTEGER(iwp) :: i, lsp
1072	INTEGER(iwp) :: ilen
1073	INTEGER(iwp) :: k
1074
1075	CHARACTER(len=16) :: spec_name
1076
1077	unit = 'illegal'
1078
1079	spec_name = TRIM(var(4:)) !< var 1:3 is 'kc_'.
1080
1081	DO lsp=1,nspec
1082	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN
1083	unit = 'ppm'
1084	ENDIF
1085	! It is possible to plant PM10 and PM25 into the gasphase chemistry code
1086	! ! as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
1087	! ! set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
1088	IF (spec_name(1:2) == 'PM') THEN
1089	unit = 'ug m-3'
1090	ENDIF
1091	ENDDO
1092
1093	DO lsp=1,nphot
1094	IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name)) THEN
1095	unit = 'sec-1'
1096	ENDIF
1097	ENDDO
1098
1099
1100	RETURN
1101	END SUBROUTINE chem_check_data_output
1102	!
1103	!------------------------------------------------------------------------------!
1104	!
1105	! Description:
1106	! ------------
1107	!> Subroutine for checking data output of profiles for chemistry model
1108	!------------------------------------------------------------------------------!
1109
1110	SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
1111
1112
1113	USE arrays_3d, &
1114	ONLY: zu
1115
1116	USE control_parameters, &
1117	ONLY: data_output_pr, message_string, air_chemistry
1118
1119	USE indices
1120
1121	USE profil_parameter
1122
1123	USE statistics
1124
1125
1126	IMPLICIT NONE
1127
1128	CHARACTER (LEN=*) :: unit !<
1129	CHARACTER (LEN=*) :: variable !<
1130	CHARACTER (LEN=*) :: dopr_unit
1131	CHARACTER(len=16) :: spec_name
1132
1133	INTEGER(iwp) :: var_count, lsp !<
1134
1135
1136	spec_name = TRIM(variable(4:))
1137	! write(9,*) 'fm #32 .. air_chemistry ', air_chemistry
1138
1139	IF ( .NOT. air_chemistry ) THEN
1140	message_string = 'data_output_pr = ' // &
1141	TRIM( data_output_pr(var_count) ) // ' is not imp' // &
1142	'lemented for air_chemistry = .FALSE.'
1143	! CALL message( 'check_parameters', 'PA0XXX', 1, 2, 0, 6, 0 )
1144	ELSE
1145	DO lsp = 1, nspec
1146	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
1147	dopr_index(var_count) = 900
1148	dopr_unit = 'ppm'
1149	hom(:,2,900,:) = SPREAD( zu, 2, statistic_regions+1 )
1150	ENDIF
1151	ENDDO
1152	ENDIF
1153
1154	END SUBROUTINE chem_check_data_output_pr
1155	!
1156	!------------------------------------------------------------------------------!
1157	!
1158	! Description:
1159	! ------------
1160	!> Subroutine defining 3D output variables for chemical species
1161	!------------------------------------------------------------------------------!
1162
1163
1164	SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1165
1166
1167	USE indices
1168
1169	USE kinds
1170
1171
1172	IMPLICIT NONE
1173
1174	CHARACTER (LEN=*) :: variable !<
1175
1176	INTEGER(iwp) :: av !<
1177	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
1178	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1179
1180	LOGICAL :: found !<
1181
1182	REAL(wp) :: fill_value !<
1183	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf
1184
1185
1186	!-- local variables
1187
1188	INTEGER :: i, j, k, lsp
1189	CHARACTER(len=16) :: spec_name
1190
1191
1192	found = .FALSE.
1193
1194	spec_name = TRIM(variable(4:))
1195	!av == 0
1196
1197	DO lsp=1,nspec
1198	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN
1199	IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM(variable) // &
1200	TRIM(chem_species(lsp)%name)
1201
1202	IF (av == 0) THEN
1203	DO i = nxl, nxr
1204	DO j = nys, nyn
1205	DO k = nzb_do, nzt_do
1206	local_pf(i,j,k) = MERGE( &
1207	chem_species(lsp)%conc(k,j,i), &
1208	REAL( fill_value, KIND = wp ), &
1209	BTEST( wall_flags_0(k,j,i), 0 ) )
1210	ENDDO
1211	ENDDO
1212	ENDDO
1213
1214	ELSE
1215	DO i = nxl, nxr
1216	DO j = nys, nyn
1217	DO k = nzb_do, nzt_do
1218	local_pf(i,j,k) = MERGE( &
1219	chem_species(lsp)%conc_av(k,j,i),&
1220	REAL( fill_value, KIND = wp ), &
1221	BTEST( wall_flags_0(k,j,i), 0 ) )
1222	ENDDO
1223	ENDDO
1224	ENDDO
1225	ENDIF
1226
1227	found = .TRUE.
1228	ENDIF
1229	ENDDO
1230
1231	RETURN
1232	END SUBROUTINE chem_data_output_3d
1233	!
1234	!------------------------------------------------------------------------------!
1235	!
1236	! Description:
1237	! ------------
1238	!> Subroutine for averaging 3D data of chemical species. Due to the fact that
1239	!> the averaged chem arrays are allocated in chem_init, no if-query concerning
1240	!> the allocation is required (in any mode). Attention: If you just specify an
1241	!> averaged output quantity in the _p3dr file during restarts the first output
1242	!> includes the time between the beginning of the restart run and the first
1243	!> output time (not necessarily the whole averaging_interval you have
1244	!> specified in your _p3d/_p3dr file )
1245	!------------------------------------------------------------------------------!
1246
1247	SUBROUTINE chem_3d_data_averaging ( mode, variable )
1248
1249	USE control_parameters
1250
1251	USE indices
1252
1253	USE kinds
1254
1255	USE surface_mod, &
1256	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
1257
1258	IMPLICIT NONE
1259
1260	CHARACTER (LEN=*) :: mode !<
1261	CHARACTER (LEN=*) :: variable !<
1262
1263	LOGICAL :: match_def !< flag indicating natural-type surface
1264	LOGICAL :: match_lsm !< flag indicating natural-type surface
1265	LOGICAL :: match_usm !< flag indicating urban-type surface
1266
1267	INTEGER(iwp) :: i !< grid index x direction
1268	INTEGER(iwp) :: j !< grid index y direction
1269	INTEGER(iwp) :: k !< grid index z direction
1270	INTEGER(iwp) :: m !< running index surface type
1271	INTEGER(iwp) :: lsp !< running index for chem spcs
1272	INTEGER(iwp) :: lsp_2 !< it looks like redundent .. will be delted ..bK
1273
1274	IF ( mode == 'allocate' ) THEN
1275	DO lsp = 1, nspec
1276	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
1277	! lsp_2 = lsp
1278	chem_species(lsp)%conc_av = 0.0_wp
1279
1280	ENDIF
1281	ENDDO
1282
1283	ELSEIF ( mode == 'sum' ) THEN
1284
1285	DO lsp = 1, nspec
1286	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
1287	! lsp_2 = lsp
1288	DO i = nxlg, nxrg
1289	DO j = nysg, nyng
1290	DO k = nzb, nzt+1
1291	chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + &
1292	chem_species(lsp)%conc(k,j,i)
1293	ENDDO
1294	ENDDO
1295	ENDDO
1296	ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') ) THEN
1297	DO i = nxl, nxr
1298	DO j = nys, nyn
1299	match_def = surf_def_h(0)%start_index(j,i) <= &
1300	surf_def_h(0)%end_index(j,i)
1301	match_lsm = surf_lsm_h%start_index(j,i) <= &
1302	surf_lsm_h%end_index(j,i)
1303	match_usm = surf_usm_h%start_index(j,i) <= &
1304	surf_usm_h%end_index(j,i)
1305
1306	IF ( match_def ) THEN
1307	m = surf_def_h(0)%end_index(j,i)
1308	chem_species(lsp)%cssws_av(j,i) = &
1309	chem_species(lsp)%cssws_av(j,i) + &
1310	surf_def_h(0)%cssws(lsp,m)
1311	ELSEIF ( match_lsm .AND. .NOT. match_usm ) THEN
1312	m = surf_lsm_h%end_index(j,i)
1313	chem_species(lsp)%cssws_av(j,i) = &
1314	chem_species(lsp)%cssws_av(j,i) + &
1315	surf_lsm_h%cssws(lsp,m)
1316	ELSEIF ( match_usm ) THEN
1317	m = surf_usm_h%end_index(j,i)
1318	chem_species(lsp)%cssws_av(j,i) = &
1319	chem_species(lsp)%cssws_av(j,i) + &
1320	surf_usm_h%cssws(lsp,m)
1321	ENDIF
1322	ENDDO
1323	ENDDO
1324	ENDIF
1325	ENDDO
1326
1327	ELSEIF ( mode == 'average' ) THEN
1328	DO lsp = 1, nspec
1329	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
1330	! lsp_2 = lsp
1331	DO i = nxlg, nxrg
1332	DO j = nysg, nyng
1333	DO k = nzb, nzt+1
1334	chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1335	ENDDO
1336	ENDDO
1337	ENDDO
1338
1339	ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') ) THEN
1340	DO i = nxlg, nxrg
1341	DO j = nysg, nyng
1342	chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
1343	ENDDO
1344	ENDDO
1345	CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )
1346	ENDIF
1347	ENDDO
1348
1349	ENDIF
1350
1351	END SUBROUTINE chem_3d_data_averaging
1352
1353	!------------------------------------------------------------------------------!
1354	!
1355	! Description:
1356	! ------------
1357	!> Subroutine to write restart data for chemistry model
1358	!------------------------------------------------------------------------------!
1359	SUBROUTINE chem_wrd_local
1360
1361
1362	IMPLICIT NONE
1363
1364	INTEGER(iwp) :: lsp !<
1365
1366	! REAL(kind=wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: chems_conc
1367
1368
1369	DO lsp = 1, nspec
1370
1371	CALL wrd_write_string( TRIM( chem_species(lsp)%name ))
1372	WRITE ( 14 ) chem_species(lsp)%conc
1373
1374	CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
1375	WRITE ( 14 ) chem_species(lsp)%conc_av
1376
1377	ENDDO
1378
1379
1380	END SUBROUTINE chem_wrd_local
1381
1382	!------------------------------------------------------------------------------!
1383	!
1384	! Description:
1385	! ------------
1386	!> Subroutine to read restart data of chemical species
1387	!------------------------------------------------------------------------------!
1388
1389	SUBROUTINE chem_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
1390	nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc, &
1391	nys_on_file, tmp_3d, found )
1392
1393	USE control_parameters
1394
1395	USE indices
1396
1397	USE pegrid
1398
1399	IMPLICIT NONE
1400
1401	CHARACTER (LEN=20) :: spc_name_av !<
1402
1403	INTEGER(iwp) :: i, lsp !<
1404	INTEGER(iwp) :: k !<
1405	INTEGER(iwp) :: nxlc !<
1406	INTEGER(iwp) :: nxlf !<
1407	INTEGER(iwp) :: nxl_on_file !<
1408	INTEGER(iwp) :: nxrc !<
1409	INTEGER(iwp) :: nxrf !<
1410	INTEGER(iwp) :: nxr_on_file !<
1411	INTEGER(iwp) :: nync !<
1412	INTEGER(iwp) :: nynf !<
1413	INTEGER(iwp) :: nyn_on_file !<
1414	INTEGER(iwp) :: nysc !<
1415	INTEGER(iwp) :: nysf !<
1416	INTEGER(iwp) :: nys_on_file !<
1417
1418	LOGICAL, INTENT(OUT) :: found
1419
1420	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
1421
1422
1423	found = .FALSE.
1424
1425
1426	IF ( ALLOCATED(chem_species) ) THEN
1427
1428	DO lsp = 1, nspec
1429
1430	!< for time-averaged chemical conc.
1431	spc_name_av = TRIM(chem_species(lsp)%name)//'_av'
1432
1433	IF (restart_string(1:length) == TRIM(chem_species(lsp)%name) ) &
1434	THEN
1435	!< read data into tmp_3d
1436	IF ( k == 1 ) READ ( 13 ) tmp_3d
1437	!< fill ..%conc in the restart run
1438	chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp, &
1439	nxlc-nbgp:nxrc+nbgp) = &
1440	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1441	found = .TRUE.
1442	ELSEIF (restart_string(1:length) == spc_name_av ) THEN
1443	IF ( k == 1 ) READ ( 13 ) tmp_3d
1444	chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp, &
1445	nxlc-nbgp:nxrc+nbgp) = &
1446	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1447	found = .TRUE.
1448	ENDIF
1449
1450	ENDDO
1451
1452	ENDIF
1453
1454
1455	END SUBROUTINE chem_rrd_local
1456
1457
1458	!------------------------------------------------------------------------------!
1459	!
1460	! Description:
1461	! ------------
1462	!> Subroutine calculating prognostic equations for chemical species
1463	!> (cache-optimized).
1464	!> Routine is called separately for each chemical species over a loop from
1465	!> prognostic_equations.
1466	!------------------------------------------------------------------------------!
1467	SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs, &
1468	i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, &
1469	flux_l_cs, diss_l_cs )
1470	USE pegrid
1471	USE advec_ws, ONLY: advec_s_ws
1472	USE advec_s_pw_mod, ONLY: advec_s_pw
1473	USE advec_s_up_mod, ONLY: advec_s_up
1474	USE diffusion_s_mod, ONLY: diffusion_s
1475	USE indices, ONLY: wall_flags_0
1476	USE surface_mod, ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
1477	surf_usm_v
1478
1479
1480	IMPLICIT NONE
1481
1482	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p, cs_scalar, tcs_scalar_m
1483
1484	INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp
1485	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: flux_s_cs !<
1486	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: diss_s_cs !<
1487	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: flux_l_cs !<
1488	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: diss_l_cs !<
1489	REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !<
1490
1491	!-- local variables
1492
1493	INTEGER :: k
1494	!
1495	!-- Tendency-terms for chem spcs.
1496	tend(:,j,i) = 0.0_wp
1497	!
1498	!-- Advection terms
1499	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1500	IF ( ws_scheme_sca ) THEN
1501	CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs, &
1502	flux_l_cs, diss_l_cs, i_omp_start, tn )
1503	ELSE
1504	CALL advec_s_pw( i, j, cs_scalar )
1505	ENDIF
1506	ELSE
1507	CALL advec_s_up( i, j, cs_scalar )
1508	ENDIF
1509
1510	!
1511
1512	!-- Diffusion terms (the last three arguments are zero)
1513
1514	CALL diffusion_s( i, j, cs_scalar, &
1515	surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:), &
1516	surf_def_h(2)%cssws(ilsp,:), &
1517	surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:), &
1518	surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:), &
1519	surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:), &
1520	surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:), &
1521	surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:), &
1522	surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:), &
1523	surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
1524
1525	!
1526	!-- Prognostic equation for chem spcs
1527	DO k = nzb+1, nzt
1528	cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d * &
1529	( tsc(2) * tend(k,j,i) + &
1530	tsc(3) * tcs_scalar_m(k,j,i) ) &
1531	- tsc(5) * rdf_sc(k) &
1532	* ( cs_scalar(k,j,i) - pr_init_cs(k) ) &
1533	) &
1534	* MERGE( 1.0_wp, 0.0_wp, &
1535	BTEST( wall_flags_0(k,j,i), 0 ) &
1536	)
1537
1538	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
1539	ENDDO
1540
1541	!
1542	!-- Calculate tendencies for the next Runge-Kutta step
1543	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1544	IF ( intermediate_timestep_count == 1 ) THEN
1545	DO k = nzb+1, nzt
1546	tcs_scalar_m(k,j,i) = tend(k,j,i)
1547	ENDDO
1548	ELSEIF ( intermediate_timestep_count < &
1549	intermediate_timestep_count_max ) THEN
1550	DO k = nzb+1, nzt
1551	tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1552	5.3125_wp * tcs_scalar_m(k,j,i)
1553	ENDDO
1554	ENDIF
1555	ENDIF
1556
1557	END SUBROUTINE chem_prognostic_equations_ij
1558
1559
1560	!------------------------------------------------------------------------------!
1561	!
1562	! Description:
1563	! ------------
1564	!> Subroutine calculating prognostic equations for chemical species
1565	!> (vector-optimized).
1566	!> Routine is called separately for each chemical species over a loop from
1567	!> prognostic_equations.
1568	!------------------------------------------------------------------------------!
1569	SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m, &
1570	pr_init_cs, ilsp )
1571
1572	USE advec_s_pw_mod, &
1573	ONLY: advec_s_pw
1574
1575	USE advec_s_up_mod, &
1576	ONLY: advec_s_up
1577
1578	USE advec_ws, &
1579	ONLY: advec_s_ws
1580
1581	USE diffusion_s_mod, &
1582	ONLY: diffusion_s
1583
1584	USE indices, &
1585	ONLY: nxl, nxr, nyn, nys, wall_flags_0
1586
1587	USE pegrid
1588
1589	USE surface_mod, &
1590	ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
1591	surf_usm_v
1592
1593	IMPLICIT NONE
1594
1595	INTEGER :: i !< running index
1596	INTEGER :: j !< running index
1597	INTEGER :: k !< running index
1598
1599	INTEGER(iwp),INTENT(IN) :: ilsp !<
1600
1601	REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !<
1602
1603	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar !<
1604	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p !<
1605	REAL(wp), DIMENSION(:,:,:), POINTER :: tcs_scalar_m !<
1606
1607
1608	!
1609	!-- Tendency terms for chemical species
1610	tend = 0.0_wp
1611	!
1612	!-- Advection terms
1613	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1614	IF ( ws_scheme_sca ) THEN
1615	CALL advec_s_ws( cs_scalar, 'kc' )
1616	ELSE
1617	CALL advec_s_pw( cs_scalar )
1618	ENDIF
1619	ELSE
1620	CALL advec_s_up( cs_scalar )
1621	ENDIF
1622	!
1623	!-- Diffusion terms (the last three arguments are zero)
1624	CALL diffusion_s( cs_scalar, &
1625	surf_def_h(0)%cssws(ilsp,:), &
1626	surf_def_h(1)%cssws(ilsp,:), &
1627	surf_def_h(2)%cssws(ilsp,:), &
1628	surf_lsm_h%cssws(ilsp,:), &
1629	surf_usm_h%cssws(ilsp,:), &
1630	surf_def_v(0)%cssws(ilsp,:), &
1631	surf_def_v(1)%cssws(ilsp,:), &
1632	surf_def_v(2)%cssws(ilsp,:), &
1633	surf_def_v(3)%cssws(ilsp,:), &
1634	surf_lsm_v(0)%cssws(ilsp,:), &
1635	surf_lsm_v(1)%cssws(ilsp,:), &
1636	surf_lsm_v(2)%cssws(ilsp,:), &
1637	surf_lsm_v(3)%cssws(ilsp,:), &
1638	surf_usm_v(0)%cssws(ilsp,:), &
1639	surf_usm_v(1)%cssws(ilsp,:), &
1640	surf_usm_v(2)%cssws(ilsp,:), &
1641	surf_usm_v(3)%cssws(ilsp,:) )
1642	!
1643	!-- Prognostic equation for chemical species
1644	DO i = nxl, nxr
1645	DO j = nys, nyn
1646	DO k = nzb+1, nzt
1647	cs_scalar_p(k,j,i) = cs_scalar(k,j,i) &
1648	+ ( dt_3d * &
1649	( tsc(2) * tend(k,j,i) &
1650	+ tsc(3) * tcs_scalar_m(k,j,i) &
1651	) &
1652	- tsc(5) * rdf_sc(k) &
1653	* ( cs_scalar(k,j,i) - pr_init_cs(k) ) &
1654	) &
1655	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1656
1657	IF ( cs_scalar_p(k,j,i) < 0.0_wp ) cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)
1658	ENDDO
1659	ENDDO
1660	ENDDO
1661	!
1662	!-- Calculate tendencies for the next Runge-Kutta step
1663	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1664	IF ( intermediate_timestep_count == 1 ) THEN
1665	DO i = nxl, nxr
1666	DO j = nys, nyn
1667	DO k = nzb+1, nzt
1668	tcs_scalar_m(k,j,i) = tend(k,j,i)
1669	ENDDO
1670	ENDDO
1671	ENDDO
1672	ELSEIF ( intermediate_timestep_count < &
1673	intermediate_timestep_count_max ) THEN
1674	DO i = nxl, nxr
1675	DO j = nys, nyn
1676	DO k = nzb+1, nzt
1677	tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i) &
1678	+ 5.3125_wp * tcs_scalar_m(k,j,i)
1679	ENDDO
1680	ENDDO
1681	ENDDO
1682	ENDIF
1683	ENDIF
1684
1685	END SUBROUTINE chem_prognostic_equations
1686
1687
1688	!------------------------------------------------------------------------------!
1689	!
1690	! Description:
1691	! ------------
1692	!> Subroutine defining header output for chemistry model
1693	!------------------------------------------------------------------------------!
1694	SUBROUTINE chem_header ( io )
1695
1696	IMPLICIT NONE
1697
1698	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
1699
1700	! print*,'the header subroutine is still not operational'
1701	!!
1702	!!-- Write chemistry model header
1703	! WRITE( io, 3 )
1704	!
1705	! IF ( radiation_scheme == "constant" ) THEN
1706	! WRITE( io, 4 ) net_radiation
1707	! ELSEIF ( radiation_scheme == "clear-sky" ) THEN
1708	! WRITE( io, 5 )
1709	! ELSEIF ( radiation_scheme == "rrtmg" ) THEN
1710	! WRITE( io, 6 )
1711	! IF ( .NOT. lw_radiation ) WRITE( io, 10 )
1712	! IF ( .NOT. sw_radiation ) WRITE( io, 11 )
1713	! ENDIF
1714	!
1715	! IF ( albedo_type == 0 ) THEN
1716	! WRITE( io, 7 ) albedo
1717	! ELSE
1718	! WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1719	! ENDIF
1720	! IF ( constant_albedo ) THEN
1721	! WRITE( io, 9 )
1722	! ENDIF
1723	!
1724	! IF ( radiation .AND. radiation_scheme /= 'constant' ) THEN
1725	! WRITE ( io, 1 ) lambda
1726	! WRITE ( io, 2 ) day_init, time_utc_init
1727	! ENDIF
1728	!
1729	! WRITE( io, 12 ) dt_radiation
1730	!
1731	! 1 FORMAT (' Geograph. longitude : lambda = ',F4.1,' degr')
1732	! 2 FORMAT (' Day of the year at model start : day_init = ',I3 &
1733	! /' UTC time at model start : time_utc_init = ',F7.1' s')
1734	! 3 FORMAT (//' Radiation model information:'/ &
1735	! ' ----------------------------'/)
1736	! 4 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, &
1737	! // 'W/m**2')
1738	! 5 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,', &
1739	! ' default)')
1740	! 6 FORMAT (' --> RRTMG scheme is used')
1741	! 7 FORMAT (/' User-specific surface albedo: albedo =', F6.3)
1742	! 8 FORMAT (/' Albedo is set for land surface type: ', A)
1743	! 9 FORMAT (/' --> Albedo is fixed during the run')
1744	!10 FORMAT (/' --> Longwave radiation is disabled')
1745	!11 FORMAT (/' --> Shortwave radiation is disabled.')
1746	!12 FORMAT (' Timestep: dt_radiation = ', F6.2, ' s')
1747	!
1748	!
1749	END SUBROUTINE chem_header
1750
1751	!------------------------------------------------------------------------------
1752	! Description:
1753	! ------------
1754	!> Subroutine reading restart data for chemistry model input parameters
1755	! (FK: To make restarts work, I had to comment this routine. We actually
1756	! don't need it, since the namelist parameters are always read in,
1757	! also in case of a restart run)
1758	!------------------------------------------------------------------------------
1759	! SUBROUTINE chem_rrd_global
1760	!
1761	! USE chem_modules
1762	!
1763	! USE control_parameters, &
1764	! ONLY: length, message_string, restart_string
1765	!
1766	!
1767	! IMPLICIT NONE
1768	!
1769	!
1770	!
1771	! DO
1772	!
1773	! SELECT CASE ( restart_string(1:length) )
1774	!
1775	! CASE ( 'bc_cs_b' )
1776	! READ ( 13 ) bc_cs_b
1777	!
1778	! CASE DEFAULT
1779	!
1780	! EXIT
1781	!
1782	! END SELECT
1783	!
1784	!!
1785	!!-- Read next string and its length
1786	! READ ( 13 ) length
1787	! READ ( 13 ) restart_string(1:length)
1788	!
1789	! ENDDO
1790	!
1791	! END SUBROUTINE chem_rrd_global
1792
1793
1794	!------------------------------------------------------------------------------!
1795	!
1796	! Description:
1797	! ------------
1798	!> Subroutine writing restart data for chemistry model input parameters
1799	! (FK: To make restarts work, I had to comment this routine. We actually
1800	! don't need it, since the namelist parameters are always read in,
1801	! also in case of a restart run)
1802	!------------------------------------------------------------------------------!
1803	! SUBROUTINE chem_wrd_global
1804	!
1805	! USE chem_modules
1806	!
1807	! USE kinds
1808	!
1809	!
1810	! IMPLICIT NONE
1811	!
1812	! INTEGER(iwp) :: lsp !< running index for chem spcs
1813	!
1814	! !
1815	! !-- Writing out input parameters that are not part of chemistry_parameters
1816	! !-- namelist (namelist parameters are anyway read in again in case of restart)
1817	! DO lsp = 1, nvar
1818	! CALL wrd_write_string( 'conc_pr_init_'//chem_species(lsp)%name )
1819	! WRITE ( 14 ) chem_species(lsp)%conc_pr_init
1820	! ENDDO
1821	!
1822	!
1823	! END SUBROUTINE chem_wrd_global
1824
1825
1826	!------------------------------------------------------------------------------!
1827	!
1828	! Description:
1829	! ------------
1830	!> Subroutine for emission
1831	!------------------------------------------------------------------------------!
1832	SUBROUTINE chem_emissions
1833
1834	USE chem_modules
1835
1836	USE netcdf_data_input_mod, &
1837	ONLY: street_type_f
1838
1839	USE surface_mod, &
1840	ONLY: surf_lsm_h
1841
1842
1843	IMPLICIT NONE
1844
1845	INTEGER(iwp) :: i !< running index for grid in x-direction
1846	INTEGER(iwp) :: j !< running index for grid in y-direction
1847	INTEGER(iwp) :: m !< running index for horizontal surfaces
1848	INTEGER(iwp) :: lsp !< running index for chem spcs
1849
1850	!
1851	!-- Comment??? (todo)
1852	IF ( street_type_f%from_file ) THEN
1853	!
1854	!-- Streets are lsm surfaces, hence, no usm surface treatment required
1855	DO m = 1, surf_lsm_h%ns
1856	i = surf_lsm_h%i(m)
1857	j = surf_lsm_h%j(m)
1858
1859	IF ( street_type_f%var(j,i) >= main_street_id .AND. &
1860	street_type_f%var(j,i) < max_street_id ) THEN
1861	DO lsp = 1, nvar
1862	surf_lsm_h%cssws(lsp,m) = emiss_factor_main * surface_csflux(lsp)
1863	ENDDO
1864	ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. &
1865	street_type_f%var(j,i) < main_street_id ) THEN
1866	DO lsp = 1, nvar
1867	surf_lsm_h%cssws(lsp,m) = emiss_factor_side * surface_csflux(lsp)
1868	ENDDO
1869	ELSE
1870	surf_lsm_h%cssws(:,m) = 0.0_wp
1871	ENDIF
1872	ENDDO
1873
1874	ENDIF
1875
1876	END SUBROUTINE chem_emissions
1877
1878
1879	END MODULE chemistry_model_mod
1880

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |