Home

Context Navigation

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

Last change on this file since 3206 was 3183, checked in by suehring, 6 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 76.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |