Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |