Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |