Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |