Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |