Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |