Home

Context Navigation

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

Last change on this file since 2912 was 2894, checked in by Giersch, 7 years ago
Reading/Writing? data in case of restart runs revised
Property svn:keywords set to `Id`
File size: 74.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |