Home

Context Navigation

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

Last change on this file since 3352 was 3318, checked in by sward, 6 years ago
Fixed faulty syntax in message string
Property svn:keywords set to `Id`
File size: 95.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 2017-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 3318 2018-10-08 11:43:01Z knoop $
29	! Fixed faulty syntax of message string
30	!
31	! 3298 2018-10-02 12:21:11Z kanani
32	! Add remarks (kanani)
33	! Merge with trunk, replaced cloud_physics by bulk_cloud_model 28.09.2018 forkel
34	! Subroutines header and chem_check_parameters added 25.09.2018 basit
35	! Removed chem_emission routine now declared in chem_emissions.f90 30.07.2018 ERUSSO
36	! Introduced emissions namelist parameters 30.07.2018 ERUSSO
37	!
38	! Timestep steering added in subroutine chem_integrate_ij and
39	! output of chosen solver in chem_parin added 30.07.2018 ketelsen
40	!
41	! chem_check_data_output_pr: added unit for PM compounds 20.07.2018 forkel
42	! replaced : by nzb+1:nzt for pt,q,ql (found by kk) 18.07.2018 forkel
43	! debugged restart run for chem species 06.07.2018 basit
44	! reorganized subroutines in alphabetical order. 27.06.2018 basit
45	! subroutine chem_parin updated for profile output 27.06.2018 basit
46	! Added humidity arrays to USE section and tmp_qvap in chem_integrate 26.6.2018 forkel
47	! Merged chemistry with with trunk (nzb_do and nzt_do in 3d output) 26.6.2018 forkel
48	!
49	! reorganized subroutines in alphabetical order. basit 22.06.2018
50	! subroutine chem_parin updated for profile output basit 22.06.2018
51	! subroutine chem_statistics added
52	! subroutine chem_check_data_output_pr add 21.06.2018 basit
53	! subroutine chem_data_output_mask added 20.05.2018 basit
54	! subroutine chem_data_output_2d added 20.05.2018 basit
55	! subroutine chem_statistics added 04.06.2018 basit
56	! subroutine chem_emissions: Set cssws to zero before setting values 20.03.2018 forkel
57	! subroutine chem_emissions: Introduced different conversion factors
58	! for PM and gaseous compounds 15.03.2018 forkel
59	! subroutine chem_emissions updated to take variable number of chem_spcs and
60	! emission factors. 13.03.2018 basit
61	! chem_boundary_conds_decycle improved. 05.03.2018 basit
62	! chem_boundary_conds_decycle subroutine added 21.02.2018 basit
63	! chem_init_profiles subroutines re-activated after correction 21.02.2018 basit
64	!
65	!
66	! 3293 2018-09-28 12:45:20Z forkel
67	! Modularization of all bulk cloud physics code components
68	!
69	! 3248 2018-09-14 09:42:06Z sward
70	! Minor formating changes
71	!
72	! 3246 2018-09-13 15:14:50Z sward
73	! Added error handling for input namelist via parin_fail_message
74	!
75	! 3241 2018-09-12 15:02:00Z raasch
76	! +nest_chemistry
77	!
78	! 3209 2018-08-27 16:58:37Z suehring
79	! Rename flags indicating outflow boundary conditions
80	!
81	! 3182 2018-07-27 13:36:03Z suehring
82	! Revise output of surface quantities in case of overhanging structures
83	!
84	! 3045 2018-05-28 07:55:41Z Giersch
85	! error messages revised
86	!
87	! 3014 2018-05-09 08:42:38Z maronga
88	! Bugfix: nzb_do and nzt_do were not used for 3d data output
89	!
90	! 3004 2018-04-27 12:33:25Z Giersch
91	! Comment concerning averaged data output added
92	!
93	! 2932 2018-03-26 09:39:22Z maronga
94	! renamed chemistry_par to chemistry_parameters
95	!
96	! 2894 2018-03-15 09:17:58Z Giersch
97	! Calculations of the index range of the subdomain on file which overlaps with
98	! the current subdomain are already done in read_restart_data_mod,
99	! chem_last_actions was renamed to chem_wrd_local, chem_read_restart_data was
100	! renamed to chem_rrd_local, chem_write_var_list was renamed to
101	! chem_wrd_global, chem_read_var_list was renamed to chem_rrd_global,
102	! chem_skip_var_list has been removed, variable named found has been
103	! introduced for checking if restart data was found, reading of restart strings
104	! has been moved completely to read_restart_data_mod, chem_rrd_local is already
105	! inside the overlap loop programmed in read_restart_data_mod, todo list has
106	! bees extended, redundant characters in chem_wrd_local have been removed,
107	! the marker * end chemistry * is not necessary anymore, strings and their
108	! respective lengths are written out and read now in case of restart runs to
109	! get rid of prescribed character lengths
110	!
111	! 2815 2018-02-19 11:29:57Z suehring
112	! Bugfix in restart mechanism,
113	! rename chem_tendency to chem_prognostic_equations,
114	! implement vector-optimized version of chem_prognostic_equations,
115	! some clean up (incl. todo list)
116	!
117	! 2773 2018-01-30 14:12:54Z suehring
118	! Declare variables required for nesting as public
119	!
120	! 2772 2018-01-29 13:10:35Z suehring
121	! Bugfix in string handling
122	!
123	! 2768 2018-01-24 15:38:29Z kanani
124	! Shorten lines to maximum length of 132 characters
125	!
126	! 2766 2018-01-22 17:17:47Z kanani
127	! Removed preprocessor directive __chem
128	!
129	! 2756 2018-01-16 18:11:14Z suehring
130	! Fill values in 3D output introduced.
131	!
132	! 2718 2018-01-02 08:49:38Z maronga
133	! Initial revision
134	!
135	!
136	!
137	!
138	! Authors:
139	! --------
140	! @author Renate Forkel
141	! @author Farah Kanani-Suehring
142	! @author Klaus Ketelsen
143	! @author Basit Khan
144	!
145	!
146	!------------------------------------------------------------------------------!
147	! Description:
148	! ------------
149	!> Chemistry model for PALM-4U
150	!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
151	!> allowed to use the chemistry model in a precursor run and additionally
152	!> not using it in a main run
153	!> @todo Update/clean-up todo list! (FK)
154	!> @todo Set proper fill values (/= 0) for chem output arrays! (FK)
155	!> @todo Add routine chem_check_parameters, add checks for inconsistent or
156	!> unallowed parameter settings.
157	!> CALL of chem_check_parameters from check_parameters. (FK)
158	!> @todo Make routine chem_header available, CALL from header.f90
159	!> (see e.g. how it is done in routine lsm_header in
160	!> land_surface_model_mod.f90). chem_header should include all setup
161	!> info about chemistry parameter settings. (FK)
162	!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
163	!> @todo Separate boundary conditions for each chem spcs to be implemented
164	!> @todo Currently only total concentration are calculated. Resolved, parameterized
165	!> and chemistry fluxes although partially and some completely coded but
166	!> are not operational/activated in this version. bK.
167	!> @todo slight differences in passive scalar and chem spcs when chem reactions
168	!> turned off. Need to be fixed. bK
169	!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
170	!> @todo chemistry error messages
171	!> @todo Format this module according to PALM coding standard (see e.g. module
172	!> template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or
173	!> D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16)
174	!
175	!------------------------------------------------------------------------------!
176
177	MODULE chemistry_model_mod
178	USE MPI
179
180	USE kinds, ONLY: wp, iwp
181	USE indices, ONLY: nz, nzb,nzt,nysg,nyng,nxlg,nxrg, &
182	nys,nyn,nx,nxl,nxr,ny,wall_flags_0
183	USE pegrid, ONLY: myid, threads_per_task
184
185	USE bulk_cloud_model_mod, &
186	ONLY: bulk_cloud_model
187
188	USE control_parameters, ONLY: bc_lr_cyc, bc_ns_cyc, dt_3d, humidity, &
189	initializing_actions, message_string, &
190	omega, tsc, intermediate_timestep_count, &
191	intermediate_timestep_count_max, &
192	timestep_scheme, use_prescribed_profile_data, ws_scheme_sca
193	USE arrays_3d, ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu
194	USE chem_gasphase_mod, ONLY: nspec, spc_names, nkppctrl, nmaxfixsteps, &
195	t_steps, chem_gasphase_integrate, vl_dim, &
196	nvar, nreact, atol, rtol, nphot, phot_names
197	USE cpulog, ONLY: cpu_log, log_point
198
199	USE chem_modules
200
201	USE statistics
202
203	IMPLICIT NONE
204	PRIVATE
205	SAVE
206
207	LOGICAL :: nest_chemistry = .TRUE. !< flag for nesting mode of chemical species, independent on parent or not
208	!
209	!- Define chemical variables
210	TYPE species_def
211	CHARACTER(LEN=8) :: name
212	CHARACTER(LEN=16) :: unit
213	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc
214	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_av
215	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_p
216	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: tconc_m
217	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: cssws_av
218	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: flux_s_cs
219	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: diss_s_cs
220	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: flux_l_cs
221	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: diss_l_cs
222	REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: conc_pr_init
223	END TYPE species_def
224
225
226	TYPE photols_def
227	CHARACTER(LEN=8) :: name
228	CHARACTER(LEN=16) :: unit
229	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: freq
230	END TYPE photols_def
231
232
233	PUBLIC species_def
234	PUBLIC photols_def
235
236
237	TYPE(species_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: chem_species
238	TYPE(photols_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: phot_frequen
239
240	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_1
241	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_2
242	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_3
243	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_av
244	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: freq_1
245
246	INTEGER,DIMENSION(nkppctrl) :: icntrl ! Fine tuning kpp
247	REAL(kind=wp),DIMENSION(nkppctrl) :: rcntrl ! Fine tuning kpp
248	LOGICAL :: decycle_chem_lr = .FALSE.
249	LOGICAL :: decycle_chem_ns = .FALSE.
250	CHARACTER (LEN=20), DIMENSION(4) :: decycle_method = &
251	(/'dirichlet','dirichlet','dirichlet','dirichlet'/)
252	!< Decycling method at horizontal boundaries,
253	!< 1=left, 2=right, 3=south, 4=north
254	!< dirichlet = initial size distribution and
255	!< chemical composition set for the ghost and
256	!< first three layers
257	!< neumann = zero gradient
258
259	REAL(kind=wp), PUBLIC :: cs_time_step
260	CHARACTER(10), PUBLIC :: photolysis_scheme
261	! 'constant',
262	! 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180
263	! 'fastj' (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
264
265	PUBLIC nest_chemistry
266	PUBLIC nspec
267	PUBLIC nreact
268	PUBLIC nvar
269	PUBLIC spc_names
270	PUBLIC spec_conc_2
271
272	!- Interface section
273	INTERFACE chem_3d_data_averaging
274	MODULE PROCEDURE chem_3d_data_averaging
275	END INTERFACE chem_3d_data_averaging
276
277	INTERFACE chem_boundary_conds
278	MODULE PROCEDURE chem_boundary_conds
279	MODULE PROCEDURE chem_boundary_conds_decycle
280	END INTERFACE chem_boundary_conds
281
282	INTERFACE chem_check_data_output
283	MODULE PROCEDURE chem_check_data_output
284	END INTERFACE chem_check_data_output
285
286	INTERFACE chem_data_output_2d
287	MODULE PROCEDURE chem_data_output_2d
288	END INTERFACE chem_data_output_2d
289
290	INTERFACE chem_data_output_3d
291	MODULE PROCEDURE chem_data_output_3d
292	END INTERFACE chem_data_output_3d
293
294	INTERFACE chem_data_output_mask
295	MODULE PROCEDURE chem_data_output_mask
296	END INTERFACE chem_data_output_mask
297
298	INTERFACE chem_check_data_output_pr
299	MODULE PROCEDURE chem_check_data_output_pr
300	END INTERFACE chem_check_data_output_pr
301
302	INTERFACE chem_check_parameters
303	MODULE PROCEDURE chem_check_parameters
304	END INTERFACE chem_check_parameters
305
306	INTERFACE chem_define_netcdf_grid
307	MODULE PROCEDURE chem_define_netcdf_grid
308	END INTERFACE chem_define_netcdf_grid
309
310	INTERFACE chem_header
311	MODULE PROCEDURE chem_header
312	END INTERFACE chem_header
313
314	INTERFACE chem_init
315	MODULE PROCEDURE chem_init
316	END INTERFACE chem_init
317
318	INTERFACE chem_init_profiles
319	MODULE PROCEDURE chem_init_profiles
320	END INTERFACE chem_init_profiles
321
322	INTERFACE chem_integrate
323	MODULE PROCEDURE chem_integrate_ij
324	END INTERFACE chem_integrate
325
326	INTERFACE chem_parin
327	MODULE PROCEDURE chem_parin
328	END INTERFACE chem_parin
329
330	INTERFACE chem_prognostic_equations
331	MODULE PROCEDURE chem_prognostic_equations
332	MODULE PROCEDURE chem_prognostic_equations_ij
333	END INTERFACE chem_prognostic_equations
334
335	INTERFACE chem_rrd_local
336	MODULE PROCEDURE chem_rrd_local
337	END INTERFACE chem_rrd_local
338
339	INTERFACE chem_statistics
340	MODULE PROCEDURE chem_statistics
341	END INTERFACE chem_statistics
342
343	INTERFACE chem_swap_timelevel
344	MODULE PROCEDURE chem_swap_timelevel
345	END INTERFACE chem_swap_timelevel
346
347	INTERFACE chem_wrd_local
348	MODULE PROCEDURE chem_wrd_local
349	END INTERFACE chem_wrd_local
350
351
352	PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
353	chem_check_data_output_pr, chem_check_parameters, &
354	chem_data_output_2d, chem_data_output_3d, chem_data_output_mask, &
355	chem_define_netcdf_grid, chem_header, chem_init, &
356	chem_init_profiles, chem_integrate, chem_parin, &
357	chem_prognostic_equations, chem_rrd_local, &
358	chem_statistics, chem_swap_timelevel, chem_wrd_local
359
360
361
362	CONTAINS
363
364	!
365	!------------------------------------------------------------------------------!
366	!
367	! Description:
368	! ------------
369	!> Subroutine for averaging 3D data of chemical species. Due to the fact that
370	!> the averaged chem arrays are allocated in chem_init, no if-query concerning
371	!> the allocation is required (in any mode). Attention: If you just specify an
372	!> averaged output quantity in the _p3dr file during restarts the first output
373	!> includes the time between the beginning of the restart run and the first
374	!> output time (not necessarily the whole averaging_interval you have
375	!> specified in your _p3d/_p3dr file )
376	!------------------------------------------------------------------------------!
377
378	SUBROUTINE chem_3d_data_averaging ( mode, variable )
379
380	USE control_parameters
381
382	USE indices
383
384	USE kinds
385
386	USE surface_mod, &
387	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
388
389	IMPLICIT NONE
390
391	CHARACTER (LEN=*) :: mode !<
392	CHARACTER (LEN=*) :: variable !<
393
394	LOGICAL :: match_def !< flag indicating natural-type surface
395	LOGICAL :: match_lsm !< flag indicating natural-type surface
396	LOGICAL :: match_usm !< flag indicating urban-type surface
397
398	INTEGER(iwp) :: i !< grid index x direction
399	INTEGER(iwp) :: j !< grid index y direction
400	INTEGER(iwp) :: k !< grid index z direction
401	INTEGER(iwp) :: m !< running index surface type
402	INTEGER(iwp) :: lsp !< running index for chem spcs
403
404
405	IF ( mode == 'allocate' ) THEN
406	DO lsp = 1, nspec
407	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
408	chem_species(lsp)%conc_av = 0.0_wp
409
410	ENDIF
411	ENDDO
412
413	ELSEIF ( mode == 'sum' ) THEN
414
415
416	DO lsp = 1, nspec
417	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
418	DO i = nxlg, nxrg
419	DO j = nysg, nyng
420	DO k = nzb, nzt+1
421	chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + &
422	chem_species(lsp)%conc(k,j,i)
423	ENDDO
424	ENDDO
425	ENDDO
426	ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') ) THEN
427	DO i = nxl, nxr
428	DO j = nys, nyn
429	match_def = surf_def_h(0)%start_index(j,i) <= &
430	surf_def_h(0)%end_index(j,i)
431	match_lsm = surf_lsm_h%start_index(j,i) <= &
432	surf_lsm_h%end_index(j,i)
433	match_usm = surf_usm_h%start_index(j,i) <= &
434	surf_usm_h%end_index(j,i)
435
436	IF ( match_def ) THEN
437	m = surf_def_h(0)%end_index(j,i)
438	chem_species(lsp)%cssws_av(j,i) = &
439	chem_species(lsp)%cssws_av(j,i) + &
440	surf_def_h(0)%cssws(lsp,m)
441	ELSEIF ( match_lsm .AND. .NOT. match_usm ) THEN
442	m = surf_lsm_h%end_index(j,i)
443	chem_species(lsp)%cssws_av(j,i) = &
444	chem_species(lsp)%cssws_av(j,i) + &
445	surf_lsm_h%cssws(lsp,m)
446	ELSEIF ( match_usm ) THEN
447	m = surf_usm_h%end_index(j,i)
448	chem_species(lsp)%cssws_av(j,i) = &
449	chem_species(lsp)%cssws_av(j,i) + &
450	surf_usm_h%cssws(lsp,m)
451	ENDIF
452	ENDDO
453	ENDDO
454	ENDIF
455	ENDDO
456
457	ELSEIF ( mode == 'average' ) THEN
458	DO lsp = 1, nspec
459	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
460	DO i = nxlg, nxrg
461	DO j = nysg, nyng
462	DO k = nzb, nzt+1
463	chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
464	ENDDO
465	ENDDO
466	ENDDO
467
468	ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') ) THEN
469	DO i = nxlg, nxrg
470	DO j = nysg, nyng
471	chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
472	ENDDO
473	ENDDO
474	CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )
475	ENDIF
476	ENDDO
477
478	ENDIF
479
480	END SUBROUTINE chem_3d_data_averaging
481
482	!
483	!------------------------------------------------------------------------------!
484	!
485	! Description:
486	! ------------
487	!> Subroutine to initialize and set all boundary conditions for chemical species
488	!------------------------------------------------------------------------------!
489
490	SUBROUTINE chem_boundary_conds( mode )
491
492	USE control_parameters, &
493	ONLY: air_chemistry, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
494	bc_radiation_s
495	USE indices, &
496	ONLY: nxl, nxr, nxlg, nxrg, nyng, nysg, nzt
497
498
499	USE arrays_3d, &
500	ONLY: dzu
501	USE surface_mod, &
502	ONLY: bc_h
503
504	CHARACTER (len=*), INTENT(IN) :: mode
505	INTEGER(iwp) :: i !< grid index x direction.
506	INTEGER(iwp) :: j !< grid index y direction.
507	INTEGER(iwp) :: k !< grid index z direction.
508	INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing.
509	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls.
510	INTEGER(iwp) :: m !< running index surface elements.
511	INTEGER(iwp) :: lsp !< running index for chem spcs.
512	INTEGER(iwp) :: lph !< running index for photolysis frequencies
513
514
515	SELECT CASE ( TRIM( mode ) )
516	CASE ( 'init' )
517
518	IF ( bc_cs_b == 'dirichlet' ) THEN
519	ibc_cs_b = 0
520	ELSEIF ( bc_cs_b == 'neumann' ) THEN
521	ibc_cs_b = 1
522	ELSE
523	message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'
524	CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 ) !<
525	ENDIF
526	!
527	!-- Set Integer flags and check for possible erroneous settings for top
528	!-- boundary condition.
529	IF ( bc_cs_t == 'dirichlet' ) THEN
530	ibc_cs_t = 0
531	ELSEIF ( bc_cs_t == 'neumann' ) THEN
532	ibc_cs_t = 1
533	ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN
534	ibc_cs_t = 2
535	ELSEIF ( bc_cs_t == 'nested' ) THEN
536	ibc_cs_t = 3
537	ELSE
538	message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'
539	CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
540	ENDIF
541
542
543	CASE ( 'set_bc_bottomtop' )
544	!-- Bottom boundary condtions for chemical species
545	DO lsp = 1, nspec
546	IF ( ibc_cs_b == 0 ) THEN
547	DO l = 0, 1
548	!-- Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
549	!-- the chem_species(nspec)%conc_p value at the topography top (k-1)
550	!-- is set; for downward-facing surfaces (l=1), kb=1, i.e. the
551	!-- value at the topography bottom (k+1) is set.
552
553	kb = MERGE( -1, 1, l == 0 )
554	!$OMP PARALLEL DO PRIVATE( i, j, k )
555	DO m = 1, bc_h(l)%ns
556	i = bc_h(l)%i(m)
557	j = bc_h(l)%j(m)
558	k = bc_h(l)%k(m)
559	chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i)
560	ENDDO
561	ENDDO
562
563	ELSEIF ( ibc_cs_b == 1 ) THEN
564	!-- in boundary_conds there is som extra loop over m here for passive tracer
565	DO l = 0, 1
566	kb = MERGE( -1, 1, l == 0 )
567	!$OMP PARALLEL DO PRIVATE( i, j, k )
568	DO m = 1, bc_h(l)%ns
569	i = bc_h(l)%i(m)
570	j = bc_h(l)%j(m)
571	k = bc_h(l)%k(m)
572	chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i)
573
574	ENDDO
575	ENDDO
576	ENDIF
577	ENDDO ! end lsp loop
578
579	!-- Top boundary conditions for chemical species - Should this not be done for all species?
580	IF ( ibc_cs_t == 0 ) THEN
581	DO lsp = 1, nspec
582	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)
583	ENDDO
584	ELSEIF ( ibc_cs_t == 1 ) THEN
585	DO lsp = 1, nspec
586	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
587	ENDDO
588	ELSEIF ( ibc_cs_t == 2 ) THEN
589	DO lsp = 1, nspec
590	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
591	ENDDO
592	ENDIF
593	!
594	CASE ( 'set_bc_lateral' )
595	!-- Lateral boundary conditions for chem species at inflow boundary
596	!-- are automatically set when chem_species concentration is
597	!-- initialized. The initially set value at the inflow boundary is not
598	!-- touched during time integration, hence, this boundary value remains
599	!-- at a constant value, which is the concentration that flows into the
600	!-- domain.
601	!-- Lateral boundary conditions for chem species at outflow boundary
602
603	IF ( bc_radiation_s ) THEN
604	DO lsp = 1, nspec
605	chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
606	ENDDO
607	ELSEIF ( bc_radiation_n ) THEN
608	DO lsp = 1, nspec
609	chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
610	ENDDO
611	ELSEIF ( bc_radiation_l ) THEN
612	DO lsp = 1, nspec
613	chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
614	ENDDO
615	ELSEIF ( bc_radiation_r ) THEN
616	DO lsp = 1, nspec
617	chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
618	ENDDO
619	ENDIF
620
621	END SELECT
622
623	END SUBROUTINE chem_boundary_conds
624
625	!
626	!------------------------------------------------------------------------------!
627	! Description:
628	! ------------
629	!> Boundary conditions for prognostic variables in chemistry: decycling in the
630	!> x-direction
631	!-----------------------------------------------------------------------------
632	SUBROUTINE chem_boundary_conds_decycle (cs_3d, cs_pr_init )
633	USE pegrid, &
634	ONLY: myid
635
636	IMPLICIT NONE
637	INTEGER(iwp) :: boundary !<
638	INTEGER(iwp) :: ee !<
639	INTEGER(iwp) :: copied !<
640	INTEGER(iwp) :: i !<
641	INTEGER(iwp) :: j !<
642	INTEGER(iwp) :: k !<
643	INTEGER(iwp) :: ss !<
644	REAL(wp), DIMENSION(nzb:nzt+1) :: cs_pr_init
645	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: cs_3d
646	REAL(wp) :: flag !< flag to mask topography grid points
647
648	flag = 0.0_wp
649
650
651	!-- Left and right boundaries
652	IF ( decycle_chem_lr .AND. bc_lr_cyc ) THEN
653
654	DO boundary = 1, 2
655
656	IF ( decycle_method(boundary) == 'dirichlet' ) THEN
657	!
658	!-- Initial profile is copied to ghost and first three layers
659	ss = 1
660	ee = 0
661	IF ( boundary == 1 .AND. nxl == 0 ) THEN
662	ss = nxlg
663	ee = nxl+2
664	ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN
665	ss = nxr-2
666	ee = nxrg
667	ENDIF
668
669	DO i = ss, ee
670	DO j = nysg, nyng
671	DO k = nzb+1, nzt
672	flag = MERGE( 1.0_wp, 0.0_wp, &
673	BTEST( wall_flags_0(k,j,i), 0 ) )
674	cs_3d(k,j,i) = cs_pr_init(k) * flag
675	ENDDO
676	ENDDO
677	ENDDO
678
679	ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN
680	!
681	!-- The value at the boundary is copied to the ghost layers to simulate
682	!-- an outlet with zero gradient
683	ss = 1
684	ee = 0
685	IF ( boundary == 1 .AND. nxl == 0 ) THEN
686	ss = nxlg
687	ee = nxl-1
688	copied = nxl
689	ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN
690	ss = nxr+1
691	ee = nxrg
692	copied = nxr
693	ENDIF
694
695	DO i = ss, ee
696	DO j = nysg, nyng
697	DO k = nzb+1, nzt
698	flag = MERGE( 1.0_wp, 0.0_wp, &
699	BTEST( wall_flags_0(k,j,i), 0 ) )
700	cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
701	ENDDO
702	ENDDO
703	ENDDO
704
705	ELSE
706	WRITE(message_string,*) &
707	'unknown decycling method: decycle_method (', &
708	boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
709	CALL message( 'chem_boundary_conds_decycle', 'CM0431', &
710	1, 2, 0, 6, 0 )
711	ENDIF
712	ENDDO
713	ENDIF
714
715
716	!-- South and north boundaries
717	IF ( decycle_chem_ns .AND. bc_ns_cyc ) THEN
718
719	DO boundary = 3, 4
720
721	IF ( decycle_method(boundary) == 'dirichlet' ) THEN
722	!
723	!-- Initial profile is copied to ghost and first three layers
724	ss = 1
725	ee = 0
726	IF ( boundary == 3 .AND. nys == 0 ) THEN
727	ss = nysg
728	ee = nys+2
729	ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN
730	ss = nyn-2
731	ee = nyng
732	ENDIF
733
734	DO i = nxlg, nxrg
735	DO j = ss, ee
736	DO k = nzb+1, nzt
737	flag = MERGE( 1.0_wp, 0.0_wp, &
738	BTEST( wall_flags_0(k,j,i), 0 ) )
739	cs_3d(k,j,i) = cs_pr_init(k) * flag
740	ENDDO
741	ENDDO
742	ENDDO
743
744
745	ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN
746	!
747	!-- The value at the boundary is copied to the ghost layers to simulate
748	!-- an outlet with zero gradient
749	ss = 1
750	ee = 0
751	IF ( boundary == 3 .AND. nys == 0 ) THEN
752	ss = nysg
753	ee = nys-1
754	copied = nys
755	ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN
756	ss = nyn+1
757	ee = nyng
758	copied = nyn
759	ENDIF
760
761	DO i = nxlg, nxrg
762	DO j = ss, ee
763	DO k = nzb+1, nzt
764	flag = MERGE( 1.0_wp, 0.0_wp, &
765	BTEST( wall_flags_0(k,j,i), 0 ) )
766	cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
767	ENDDO
768	ENDDO
769	ENDDO
770
771	ELSE
772	WRITE(message_string,*) &
773	'unknown decycling method: decycle_method (', &
774	boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
775	CALL message( 'chem_boundary_conds_decycle', 'CM0432', &
776	1, 2, 0, 6, 0 )
777	ENDIF
778	ENDDO
779	ENDIF
780	END SUBROUTINE chem_boundary_conds_decycle
781	!
782	!------------------------------------------------------------------------------!
783	!
784	! Description:
785	! ------------
786	!> Subroutine for checking data output for chemical species
787	!------------------------------------------------------------------------------!
788
789	SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
790
791
792	USE control_parameters, &
793	ONLY: data_output, message_string
794
795	IMPLICIT NONE
796
797	CHARACTER (LEN=*) :: unit !<
798	CHARACTER (LEN=*) :: var !<
799
800	INTEGER(iwp) :: i, lsp
801	INTEGER(iwp) :: ilen
802	INTEGER(iwp) :: k
803
804	CHARACTER(len=16) :: spec_name
805
806	unit = 'illegal'
807
808	spec_name = TRIM(var(4:)) !< var 1:3 is 'kc_'.
809
810	DO lsp=1,nspec
811	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN
812	unit = 'ppm'
813	ENDIF
814	! It is possible to plant PM10 and PM25 into the gasphase chemistry code
815	! as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
816	! set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
817	IF (spec_name(1:2) == 'PM') THEN
818	unit = 'kg m-3'
819	ENDIF
820	ENDDO
821
822	DO lsp=1,nphot
823	IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name)) THEN
824	unit = 'sec-1'
825	ENDIF
826	ENDDO
827
828
829	RETURN
830	END SUBROUTINE chem_check_data_output
831	!
832	!------------------------------------------------------------------------------!
833	!
834	! Description:
835	! ------------
836	!> Subroutine for checking data output of profiles for chemistry model
837	!------------------------------------------------------------------------------!
838
839	SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
840
841	USE arrays_3d
842	USE control_parameters, &
843	ONLY: data_output_pr, message_string, air_chemistry
844	USE indices
845	USE profil_parameter
846	USE statistics
847
848
849	IMPLICIT NONE
850
851	CHARACTER (LEN=*) :: unit !<
852	CHARACTER (LEN=*) :: variable !<
853	CHARACTER (LEN=*) :: dopr_unit
854	CHARACTER(len=16) :: spec_name
855
856	INTEGER(iwp) :: var_count, lsp !<
857
858
859	spec_name = TRIM(variable(4:))
860
861	IF ( .NOT. air_chemistry ) THEN
862	message_string = 'data_output_pr = ' // &
863	TRIM( data_output_pr(var_count) ) // ' is not imp' // &
864	'lemented for air_chemistry = .FALSE.'
865	CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )
866
867	ELSE
868	DO lsp = 1, nspec
869	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
870	cs_pr_count = cs_pr_count+1
871	cs_pr_index(cs_pr_count) = lsp
872	dopr_index(var_count) = pr_palm+cs_pr_count
873	dopr_unit = 'ppm'
874	IF (spec_name(1:2) == 'PM') THEN
875	dopr_unit = 'kg m-3'
876	ENDIF
877	hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
878	unit = dopr_unit
879	ENDIF
880	ENDDO
881	ENDIF
882
883	END SUBROUTINE chem_check_data_output_pr
884
885	!
886	!------------------------------------------------------------------------------!
887	! Description:
888	! ------------
889	!> Check parameters routine for chemistry_model_mod
890	!------------------------------------------------------------------------------!
891	SUBROUTINE chem_check_parameters
892
893	IMPLICIT NONE
894
895	LOGICAL :: found
896	INTEGER (iwp) :: lsp_usr !< running index for user defined chem spcs
897	INTEGER (iwp) :: lsp !< running index for chem spcs.
898
899
900	!!-- check for chemical reactions status
901	IF ( chem_gasphase_on ) THEN
902	message_string = 'Chemical reactions: ON'
903	CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
904	ELSEIF ( .not. (chem_gasphase_on) ) THEN
905	message_string = 'Chemical reactions: OFF'
906	CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
907	ENDIF
908
909	!-- check for chemistry time-step
910	IF ( call_chem_at_all_substeps ) THEN
911	message_string = 'Chemistry is calculated at all meteorology time-step'
912	CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
913	ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN
914	message_string = 'Sub-time-steps are skipped for chemistry time-steps'
915	CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
916	ENDIF
917
918	!-- check for photolysis scheme
919	IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant') ) THEN
920	message_string = 'Incorrect photolysis scheme selected, please check spelling'
921	CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
922	ENDIF
923
924	!-- check for decycling of chem species
925	IF ( (.not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) ) THEN
926	message_string = 'Incorrect boundary conditions. Only neumann or ' &
927	// 'dirichlet &available for decycling chemical species '
928	CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
929	ENDIF
930
931	!-- check for initial chem species input
932	lsp_usr = 1
933	lsp = 1
934	DO WHILE ( cs_name (lsp_usr) /= 'novalue')
935	found = .FALSE.
936	DO lsp = 1, nvar
937	IF ( TRIM(cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
938	found = .TRUE.
939	EXIT
940	ENDIF
941	ENDDO
942	IF ( .not. found ) THEN
943	message_string = 'Incorrect input for initial surface vaue: ' // TRIM(cs_name(lsp_usr))
944	CALL message( 'chem_check_parameters', 'CM0427', 0, 1, 0, 6, 0 )
945	ENDIF
946	lsp_usr = lsp_usr + 1
947	ENDDO
948
949	!-- check for surface emission flux chem species
950
951	lsp_usr = 1
952	lsp = 1
953	DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
954	found = .FALSE.
955	DO lsp = 1, nvar
956	IF ( TRIM(surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
957	found = .TRUE.
958	EXIT
959	ENDIF
960	ENDDO
961	IF ( .not. found ) THEN
962	message_string = 'Incorrect input of chemical species for surface emission fluxes: ' &
963	// TRIM(surface_csflux_name(lsp_usr))
964	CALL message( 'chem_check_parameters', 'CM0428', 0, 1, 0, 6, 0 )
965	ENDIF
966	lsp_usr = lsp_usr + 1
967	ENDDO
968
969	END SUBROUTINE chem_check_parameters
970
971	!
972	!------------------------------------------------------------------------------!
973	!
974	! Description:
975	! ------------
976	!> Subroutine defining 2D output variables for chemical species
977	!> @todo: Remove "mode" from argument list, not used.
978	!------------------------------------------------------------------------------!
979
980	SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf, &
981	two_d, nzb_do, nzt_do, fill_value )
982
983	USE indices
984
985	USE kinds
986
987	USE pegrid, ONLY: myid, threads_per_task
988
989	IMPLICIT NONE
990
991	CHARACTER (LEN=*) :: grid !<
992	CHARACTER (LEN=*) :: mode !<
993	CHARACTER (LEN=*) :: variable !<
994	INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data
995	INTEGER(iwp) :: nzb_do !< lower limit of the domain (usually nzb)
996	INTEGER(iwp) :: nzt_do !< upper limit of the domain (usually nzt+1)
997	LOGICAL :: found !<
998	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
999	REAL(wp) :: fill_value
1000	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf !<
1001
1002	!-- local variables.
1003	CHARACTER(len=16) :: spec_name
1004	INTEGER(iwp) :: lsp
1005	INTEGER(iwp) :: i !< grid index along x-direction
1006	INTEGER(iwp) :: j !< grid index along y-direction
1007	INTEGER(iwp) :: k !< grid index along z-direction
1008	INTEGER(iwp) :: m !< running index surface elements
1009	INTEGER(iwp) :: char_len !< length of a character string
1010	found = .TRUE.
1011	char_len = LEN_TRIM(variable)
1012
1013	spec_name = TRIM( variable(4:char_len-3) )
1014
1015	DO lsp=1,nspec
1016	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) .AND. &
1017	( (variable(char_len-2:) == '_xy') .OR. &
1018	(variable(char_len-2:) == '_xz') .OR. &
1019	(variable(char_len-2:) == '_yz') ) ) THEN
1020
1021	IF(myid == 0) WRITE(6,*) 'Output of species ' // TRIM(variable) // &
1022	TRIM(chem_species(lsp)%name)
1023	IF (av == 0) THEN
1024	DO i = nxl, nxr
1025	DO j = nys, nyn
1026	DO k = nzb_do, nzt_do
1027	local_pf(i,j,k) = MERGE( &
1028	chem_species(lsp)%conc(k,j,i), &
1029	REAL( fill_value, KIND = wp ), &
1030	BTEST( wall_flags_0(k,j,i), 0 ) )
1031
1032
1033	ENDDO
1034	ENDDO
1035	ENDDO
1036
1037	ELSE
1038	DO i = nxl, nxr
1039	DO j = nys, nyn
1040	DO k = nzb_do, nzt_do
1041	local_pf(i,j,k) = MERGE( &
1042	chem_species(lsp)%conc(k,j,i), &
1043	REAL( fill_value, KIND = wp ), &
1044	BTEST( wall_flags_0(k,j,i), 0 ) )
1045	ENDDO
1046	ENDDO
1047	ENDDO
1048	ENDIF
1049	grid = 'zu'
1050	ENDIF
1051	ENDDO
1052
1053	RETURN
1054
1055	END SUBROUTINE chem_data_output_2d
1056
1057	!
1058	!------------------------------------------------------------------------------!
1059	!
1060	! Description:
1061	! ------------
1062	!> Subroutine defining 3D output variables for chemical species
1063	!------------------------------------------------------------------------------!
1064
1065	SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1066
1067
1068	USE indices
1069
1070	USE kinds
1071
1072
1073	IMPLICIT NONE
1074
1075	CHARACTER (LEN=*) :: variable !<
1076	INTEGER(iwp) :: av !<
1077	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
1078	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1079
1080	LOGICAL :: found !<
1081
1082	REAL(wp) :: fill_value !<
1083	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf
1084
1085
1086	!-- local variables
1087
1088	INTEGER :: i, j, k, lsp
1089	CHARACTER(len=16) :: spec_name
1090
1091
1092	found = .FALSE.
1093
1094	spec_name = TRIM(variable(4:))
1095
1096	DO lsp=1,nspec
1097	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN
1098	IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM(variable) // &
1099	TRIM(chem_species(lsp)%name)
1100
1101	IF (av == 0) THEN
1102	DO i = nxl, nxr
1103	DO j = nys, nyn
1104	DO k = nzb_do, nzt_do
1105	local_pf(i,j,k) = MERGE( &
1106	chem_species(lsp)%conc(k,j,i), &
1107	REAL( fill_value, KIND = wp ), &
1108	BTEST( wall_flags_0(k,j,i), 0 ) )
1109	ENDDO
1110	ENDDO
1111	ENDDO
1112
1113	ELSE
1114	DO i = nxl, nxr
1115	DO j = nys, nyn
1116	DO k = nzb_do, nzt_do
1117	local_pf(i,j,k) = MERGE( &
1118	chem_species(lsp)%conc_av(k,j,i),&
1119	REAL( fill_value, KIND = wp ), &
1120	BTEST( wall_flags_0(k,j,i), 0 ) )
1121	ENDDO
1122	ENDDO
1123	ENDDO
1124	ENDIF
1125
1126	found = .TRUE.
1127	ENDIF
1128	ENDDO
1129
1130	RETURN
1131	END SUBROUTINE chem_data_output_3d
1132	!
1133	!------------------------------------------------------------------------------!
1134	!
1135	! Description:
1136	! ------------
1137	!> Subroutine defining mask output variables for chemical species
1138	!------------------------------------------------------------------------------!
1139
1140	SUBROUTINE chem_data_output_mask( av, variable, found, local_pf )
1141
1142	USE control_parameters
1143	USE indices
1144	USE kinds
1145	USE pegrid, ONLY: myid, threads_per_task
1146
1147
1148	IMPLICIT NONE
1149
1150	CHARACTER (LEN=*):: variable !<
1151	INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data
1152	LOGICAL :: found !<
1153	REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
1154	local_pf !<
1155
1156
1157	!-- local variables.
1158	CHARACTER(len=16) :: spec_name
1159	INTEGER(iwp) :: lsp
1160	INTEGER(iwp) :: i !< grid index along x-direction
1161	INTEGER(iwp) :: j !< grid index along y-direction
1162	INTEGER(iwp) :: k !< grid index along z-direction
1163	found = .TRUE.
1164
1165	spec_name = TRIM( variable(4:) )
1166	!av == 0
1167
1168	DO lsp=1,nspec
1169	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) ) THEN
1170
1171	IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM(variable) // &
1172	TRIM(chem_species(lsp)%name)
1173	IF (av == 0) THEN
1174	DO i = 1, mask_size_l(mid,1)
1175	DO j = 1, mask_size_l(mid,2)
1176	DO k = 1, mask_size(mid,3)
1177	local_pf(i,j,k) = chem_species(lsp)%conc(mask_k(mid,k), &
1178	mask_j(mid,j), mask_i(mid,i))
1179	ENDDO
1180	ENDDO
1181	ENDDO
1182
1183	ELSE
1184	DO i = 1, mask_size_l(mid,1)
1185	DO j = 1, mask_size_l(mid,2)
1186	DO k = 1, mask_size_l(mid,3)
1187	local_pf(i,j,k) = chem_species(lsp)%conc_av(mask_k(mid,k), &
1188	mask_j(mid,j), mask_i(mid,i))
1189	ENDDO
1190	ENDDO
1191	ENDDO
1192	ENDIF
1193	found = .FALSE.
1194	ENDIF
1195	ENDDO
1196
1197	RETURN
1198
1199	END SUBROUTINE chem_data_output_mask
1200
1201	!
1202	!------------------------------------------------------------------------------!
1203	!
1204	! Description:
1205	! ------------
1206	!> Subroutine defining appropriate grid for netcdf variables.
1207	!> It is called out from subroutine netcdf.
1208	!------------------------------------------------------------------------------!
1209	SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1210
1211	IMPLICIT NONE
1212
1213	CHARACTER (LEN=*), INTENT(IN) :: var !<
1214	LOGICAL, INTENT(OUT) :: found !<
1215	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
1216	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
1217	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
1218
1219	found = .TRUE.
1220
1221	IF (var(1:3) == 'kc_') THEN !< always the same grid for chemistry variables
1222	grid_x = 'x'
1223	grid_y = 'y'
1224	grid_z = 'zu'
1225	ELSE
1226	found = .FALSE.
1227	grid_x = 'none'
1228	grid_y = 'none'
1229	grid_z = 'none'
1230	ENDIF
1231
1232
1233	END SUBROUTINE chem_define_netcdf_grid
1234	!
1235	!------------------------------------------------------------------------------!
1236	!
1237	! Description:
1238	! ------------
1239	!> Subroutine defining header output for chemistry model
1240	!------------------------------------------------------------------------------!
1241	SUBROUTINE chem_header ( io )
1242
1243	IMPLICIT NONE
1244
1245	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
1246	INTEGER(iwp) :: lsp !< running index for chem spcs
1247	INTEGER(iwp) :: cs_fixed
1248	CHARACTER (LEN=80) :: docsflux_chr
1249	CHARACTER (LEN=80) :: docsinit_chr
1250
1251	!
1252	!-- Write chemistry model header
1253	WRITE( io, 1 )
1254
1255	!-- Gasphase reaction status
1256	IF ( chem_gasphase_on ) THEN
1257	WRITE( io, 2 )
1258	ELSE
1259	WRITE( io, 3 )
1260	ENDIF
1261
1262	! Chemistry time-step
1263	WRITE ( io, 4 ) cs_time_step
1264
1265	!-- Emission mode info
1266	IF ( mode_emis == "DEFAULT" ) THEN
1267	WRITE( io, 5 )
1268	ELSEIF ( mode_emis == "PARAMETERIZED" ) THEN
1269	WRITE( io, 6 )
1270	ELSEIF ( mode_emis == "PRE-PROCESSED" ) THEN
1271	WRITE( io, 7 )
1272	ENDIF
1273
1274	!-- Photolysis scheme info
1275	IF ( photolysis_scheme == "simple" ) THEN
1276	WRITE( io, 8 )
1277	ELSEIF (photolysis_scheme == "conastant" ) THEN
1278	WRITE( io, 9 )
1279	ENDIF
1280
1281	!-- Emission flux info
1282	lsp = 1
1283	docsflux_chr ='Chemical species for surface emission flux: '
1284	DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
1285	docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ','
1286	IF ( LEN_TRIM( docsflux_chr ) >= 75 ) THEN
1287	WRITE ( io, 10 ) docsflux_chr
1288	docsflux_chr = ' '
1289	ENDIF
1290	lsp = lsp + 1
1291	ENDDO
1292
1293	IF ( docsflux_chr /= '' ) THEN
1294	WRITE ( io, 10 ) docsflux_chr
1295	ENDIF
1296
1297
1298	!-- initializatoin of Surface and profile chemical species
1299
1300	lsp = 1
1301	docsinit_chr ='Chemical species for initial surface and profile emissions: '
1302	DO WHILE ( cs_name(lsp) /= 'novalue' )
1303	docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ','
1304	IF ( LEN_TRIM( docsinit_chr ) >= 75 ) THEN
1305	WRITE ( io, 11 ) docsinit_chr
1306	docsinit_chr = ' '
1307	ENDIF
1308	lsp = lsp + 1
1309	ENDDO
1310
1311	IF ( docsinit_chr /= '' ) THEN
1312	WRITE ( io, 11 ) docsinit_chr
1313	ENDIF
1314
1315	!-- number of variable and fix chemical species and number of reactions
1316	cs_fixed = nspec - nvar
1317	WRITE ( io, * ) ' --> Chemical species, variable: ', nvar
1318	WRITE ( io, * ) ' --> Chemical species, fixed : ', cs_fixed
1319	WRITE ( io, * ) ' --> Total number of reactions : ', nreact
1320
1321
1322	1 FORMAT (//' Chemistry model information:'/ &
1323	' ----------------------------'/)
1324	2 FORMAT (' --> Chemical reactions are turned on')
1325	3 FORMAT (' --> Chemical reactions are turned off')
1326	4 FORMAT (' --> Time-step for chemical species: ',F6.2, ' s')
1327	5 FORMAT (' --> Emission mode = DEFAULT ')
1328	6 FORMAT (' --> Emission mode = PARAMETERIZED ')
1329	7 FORMAT (' --> Emission mode = PRE-PROCESSED ')
1330	8 FORMAT (' --> Photolysis scheme used = simple ')
1331	9 FORMAT (' --> Photolysis scheme used = constant ')
1332	10 FORMAT (/' ',A)
1333	11 FORMAT (/' ',A)
1334	!
1335	!
1336	END SUBROUTINE chem_header
1337
1338	!
1339	!------------------------------------------------------------------------------!
1340	!
1341	! Description:
1342	! ------------
1343	!> Subroutine initializating chemistry_model_mod
1344	!------------------------------------------------------------------------------!
1345	SUBROUTINE chem_init
1346
1347	USE control_parameters, &
1348	ONLY: message_string, io_blocks, io_group, turbulent_inflow
1349	USE arrays_3d, &
1350	ONLY: mean_inflow_profiles
1351	USE pegrid
1352
1353	IMPLICIT none
1354	!-- local variables
1355	INTEGER :: i,j !< running index for for horiz numerical grid points
1356	INTEGER :: lsp !< running index for chem spcs
1357	INTEGER :: lpr_lev !< running index for chem spcs profile level
1358	!
1359	!-- NOPOINTER version not implemented yet
1360	! #if defined( __nopointer )
1361	! message_string = 'The chemistry module only runs with POINTER version'
1362	! CALL message( 'chemistry_model_mod', 'CM0001', 1, 2, 0, 6, 0 )
1363	! #endif
1364	!
1365	!-- Allocate memory for chemical species
1366	ALLOCATE( chem_species(nspec) )
1367	ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1368	ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1369	ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1370	ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1371	ALLOCATE( phot_frequen(nphot) )
1372	ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
1373	ALLOCATE( bc_cs_t_val(nspec) )
1374	!
1375	!-- Initialize arrays
1376	spec_conc_1 (:,:,:,:) = 0.0_wp
1377	spec_conc_2 (:,:,:,:) = 0.0_wp
1378	spec_conc_3 (:,:,:,:) = 0.0_wp
1379	spec_conc_av(:,:,:,:) = 0.0_wp
1380
1381
1382	DO lsp = 1, nspec
1383	chem_species(lsp)%name = spc_names(lsp)
1384
1385	chem_species(lsp)%conc (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1 (:,:,:,lsp)
1386	chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2 (:,:,:,lsp)
1387	chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_3 (:,:,:,lsp)
1388	chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_av(:,:,:,lsp)
1389
1390	ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))
1391	chem_species(lsp)%cssws_av = 0.0_wp
1392	!
1393	!-- The following block can be useful when emission module is not applied. &
1394	!-- if emission module is applied the following block will be overwritten.
1395	ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))
1396	ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))
1397	ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
1398	ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
1399	chem_species(lsp)%flux_s_cs = 0.0_wp
1400	chem_species(lsp)%flux_l_cs = 0.0_wp
1401	chem_species(lsp)%diss_s_cs = 0.0_wp
1402	chem_species(lsp)%diss_l_cs = 0.0_wp
1403	!
1404	!-- Allocate memory for initial concentration profiles
1405	!-- (concentration values come from namelist)
1406	!-- (@todo (FK): Because of this, chem_init is called in palm before
1407	!-- check_parameters, since conc_pr_init is used there.
1408	!-- We have to find another solution since chem_init should
1409	!-- eventually be called from init_3d_model!!)
1410	ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
1411	chem_species(lsp)%conc_pr_init(:) = 0.0_wp
1412
1413	ENDDO
1414
1415	!
1416	!-- Initial concentration of profiles is prescribed by parameters cs_profile
1417	!-- and cs_heights in the namelist &chemistry_parameters
1418	CALL chem_init_profiles
1419
1420
1421	!
1422	!-- Initialize model variables
1423
1424
1425	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1426	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
1427
1428
1429	!-- First model run of a possible job queue.
1430	!-- Initial profiles of the variables must be computed.
1431	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
1432	CALL location_message( 'initializing with 1D chemistry model profiles', .FALSE. )
1433	!
1434	!-- Transfer initial profiles to the arrays of the 3D model
1435	DO lsp = 1, nspec
1436	DO i = nxlg, nxrg
1437	DO j = nysg, nyng
1438	DO lpr_lev = 1, nz + 1
1439	chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
1440	ENDDO
1441	ENDDO
1442	ENDDO
1443	ENDDO
1444
1445	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
1446	THEN
1447	CALL location_message( 'initializing with constant chemistry profiles', .FALSE. )
1448
1449	DO lsp = 1, nspec
1450	DO i = nxlg, nxrg
1451	DO j = nysg, nyng
1452	chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init
1453	ENDDO
1454	ENDDO
1455	ENDDO
1456
1457	ENDIF
1458
1459	!
1460	!-- If required, change the surface chem spcs at the start of the 3D run
1461	IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN
1462	DO lsp = 1, nspec
1463	chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) + &
1464	cs_surface_initial_change(lsp)
1465	ENDDO
1466	ENDIF
1467	!
1468	!-- Initiale old and new time levels.
1469	DO lsp = 1, nvar
1470	chem_species(lsp)%tconc_m = 0.0_wp
1471	chem_species(lsp)%conc_p = chem_species(lsp)%conc
1472	ENDDO
1473
1474	ENDIF
1475
1476
1477
1478	!--- new code add above this line
1479	DO lsp = 1, nphot
1480	phot_frequen(lsp)%name = phot_names(lsp)
1481	! IF( myid == 0 ) THEN
1482	! WRITE(6,'(a,i4,3x,a)') 'Photolysis: ',lsp,trim(phot_names(lsp))
1483	! ENDIF
1484	phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => freq_1(:,:,:,lsp)
1485	ENDDO
1486
1487	RETURN
1488
1489	END SUBROUTINE chem_init
1490
1491	!
1492	!------------------------------------------------------------------------------!
1493	!
1494	! Description:
1495	! ------------
1496	!> Subroutine defining initial vertical profiles of chemical species (given by
1497	!> namelist parameters chem_profiles and chem_heights) --> which should work
1498	!> analogue to parameters u_profile, v_profile and uv_heights)
1499	!------------------------------------------------------------------------------!
1500	SUBROUTINE chem_init_profiles !< SUBROUTINE is called from chem_init in case of
1501	!< TRIM( initializing_actions ) /= 'read_restart_data'
1502	!< We still need to see what has to be done in case of restart run
1503	USE chem_modules
1504
1505	IMPLICIT NONE
1506
1507	!-- Local variables
1508	INTEGER :: lsp !< running index for number of species in derived data type species_def
1509	INTEGER :: lsp_usr !< running index for number of species (user defined) in cs_names, cs_profiles etc
1510	INTEGER :: lpr_lev !< running index for profile level for each chem spcs.
1511	INTEGER :: npr_lev !< the next available profile lev
1512
1513	!-----------------
1514	!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
1515	!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
1516	!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
1517	!-- "cs_heights" are prescribed, their values will!override the constant profile given by
1518	!-- "cs_surface".
1519	!
1520	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
1521	lsp_usr = 1
1522	DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) !'novalue' is the default
1523	DO lsp = 1, nspec !
1524	!-- create initial profile (conc_pr_init) for each chemical species
1525	IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) ) THEN !
1526	IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN
1527	!-- set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
1528	DO lpr_lev = 0, nzt+1
1529	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
1530	ENDDO
1531	ELSE
1532	IF ( cs_heights(1,1) /= 0.0_wp ) THEN
1533	message_string = 'The surface value of cs_heights must be 0.0'
1534	CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
1535	ENDIF
1536
1537	use_prescribed_profile_data = .TRUE.
1538
1539	npr_lev = 1
1540	! chem_species(lsp)%conc_pr_init(0) = 0.0_wp
1541	DO lpr_lev = 1, nz+1
1542	IF ( npr_lev < 100 ) THEN
1543	DO WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
1544	npr_lev = npr_lev + 1
1545	IF ( npr_lev == 100 ) THEN
1546	message_string = 'number of chem spcs exceeding the limit'
1547	CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )
1548	EXIT
1549	ENDIF
1550	ENDDO
1551	ENDIF
1552	IF ( npr_lev < 100 .AND. cs_heights(lsp_usr, npr_lev + 1) /= 9999999.9_wp ) THEN
1553	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) + &
1554	( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) / &
1555	( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) * &
1556	( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
1557	ELSE
1558	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
1559	ENDIF
1560	ENDDO
1561	ENDIF
1562	!-- If a profile is prescribed explicity using cs_profiles and cs_heights, then
1563	!-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
1564	!-- on the cs_profiles(lsp_usr,:) and cs_heights(lsp_usr,:).
1565	ENDIF
1566	ENDDO
1567	lsp_usr = lsp_usr + 1
1568	ENDDO
1569	ENDIF
1570
1571	END SUBROUTINE chem_init_profiles
1572	!
1573	!------------------------------------------------------------------------------!
1574	!
1575	! Description:
1576	! ------------
1577	!> Subroutine to integrate chemical species in the given chemical mechanism
1578	!------------------------------------------------------------------------------!
1579
1580	SUBROUTINE chem_integrate_ij (i, j)
1581
1582	USE cpulog, &
1583	ONLY: cpu_log, log_point
1584	USE statistics, &
1585	ONLY: weight_pres
1586	USE control_parameters, &
1587	ONLY: dt_3d, intermediate_timestep_count,simulated_time
1588
1589	IMPLICIT none
1590	INTEGER,INTENT(IN) :: i,j
1591
1592	!-- local variables
1593	INTEGER(iwp) :: lsp !< running index for chem spcs.
1594	INTEGER(iwp) :: lph !< running index for photolysis frequencies
1595	INTEGER :: k,m,istatf
1596	INTEGER,DIMENSION(20) :: istatus
1597	REAL(kind=wp),DIMENSION(nzb+1:nzt,nspec) :: tmp_conc
1598	REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_temp
1599	REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_qvap
1600	REAL(kind=wp),DIMENSION(nzb+1:nzt,nphot) :: tmp_phot
1601	REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_fact
1602	REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_fact_i !< conversion factor between
1603	!< molecules cm^{-3} and ppm
1604
1605	INTEGER,DIMENSION(nzb+1:nzt) :: nacc !< Number of accepted steps
1606	INTEGER,DIMENSION(nzb+1:nzt) :: nrej !< Number of rejected steps
1607
1608	REAL(wp) :: conv !< conversion factor
1609	REAL(wp), PARAMETER :: ppm2fr = 1.0e-6_wp !< Conversion factor ppm to fraction
1610	REAL(wp), PARAMETER :: xm_air = 28.96_wp !< Mole mass of dry air
1611	REAL(wp), PARAMETER :: xm_h2o = 18.01528_wp !< Mole mass of water vapor
1612	REAL(wp), PARAMETER :: pref_i = 1._wp / 100000.0_wp !< inverse reference pressure (1/Pa)
1613	REAL(wp), PARAMETER :: t_std = 273.15_wp !< standard pressure (Pa)
1614	REAL(wp), PARAMETER :: p_std = 101325.0_wp !< standard pressure (Pa)
1615	REAL(wp), PARAMETER :: vmolcm = 22.414e3_wp !< Mole volume (22.414 l) in cm^{-3}
1616	REAL(wp), PARAMETER :: xna = 6.022e23_wp !< Avogadro number (molecules/mol)
1617
1618	REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
1619
1620
1621	REAL(kind=wp) :: dt_chem
1622
1623	CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' )
1624	!< set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
1625	IF (chem_gasphase_on) THEN
1626	nacc = 0
1627	nrej = 0
1628
1629	tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
1630	! ppm to molecules/cm**3
1631	! tmp_fact = 1.e-6_wp6.022e23_wp/(22.414_wp1000._wp) * 273.15_wp * hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp )
1632	conv = ppm2fr * xna / vmolcm
1633	tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
1634	tmp_fact_i = 1.0_wp/tmp_fact
1635
1636	IF ( humidity ) THEN
1637	IF ( bulk_cloud_model ) THEN
1638	tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) * xm_air/xm_h2o * tmp_fact(:)
1639	ELSE
1640	tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * tmp_fact(:)
1641	ENDIF
1642	ELSE
1643	tmp_qvap(:) = 0.01 * tmp_fact(:) !< Constant value for q if water vapor is not computed
1644	ENDIF
1645
1646	DO lsp = 1,nspec
1647	tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
1648	ENDDO
1649
1650	DO lph = 1,nphot
1651	tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)
1652	ENDDO
1653
1654	IF(myid == 0 .AND. chem_debug0 ) THEN
1655	IF (i == 10 .and. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
1656	ENDIF
1657
1658	!-- Compute length of time step
1659	IF ( call_chem_at_all_substeps ) THEN
1660	dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
1661	ELSE
1662	dt_chem = dt_3d
1663	ENDIF
1664
1665	cs_time_step = dt_chem
1666
1667	CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' )
1668
1669	IF(maxval(rcntrl) > 0.0) THEN ! Only if rcntrl is set
1670	IF( simulated_time <= 2*dt_3d) THEN
1671	rcntrl_local = 0
1672	WRITE(9,'(a,2f10.3)') 'USE Default rcntrl in the first steps ',simulated_time,dt_3d
1673	ELSE
1674	rcntrl_local = rcntrl
1675	ENDIF
1676	ELSE
1677	rcntrl_local = 0
1678	END IF
1679
1680	CALL chem_gasphase_integrate (dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
1681	icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus)
1682
1683	CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' )
1684
1685	DO lsp = 1,nspec
1686	chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
1687	ENDDO
1688
1689
1690	ENDIF
1691	CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' )
1692
1693	RETURN
1694	END SUBROUTINE chem_integrate_ij
1695	!
1696	!------------------------------------------------------------------------------!
1697	!
1698	! Description:
1699	! ------------
1700	!> Subroutine defining parin for &chemistry_parameters for chemistry model
1701	!------------------------------------------------------------------------------!
1702	SUBROUTINE chem_parin
1703
1704	USE chem_modules
1705	USE control_parameters
1706
1707	USE kinds
1708	USE pegrid
1709	USE statistics
1710
1711	IMPLICIT NONE
1712
1713	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
1714	CHARACTER (LEN=3) :: cs_prefix
1715
1716	REAL(wp), DIMENSION(nmaxfixsteps) :: my_steps !< List of fixed timesteps my_step(1) = 0.0 automatic stepping
1717	INTEGER(iwp) :: i !<
1718	INTEGER(iwp) :: j !<
1719	INTEGER(iwp) :: max_pr_cs_tmp !<
1720
1721
1722	NAMELIST /chemistry_parameters/ bc_cs_b, &
1723	bc_cs_t, &
1724	call_chem_at_all_substeps, &
1725	chem_debug0, &
1726	chem_debug1, &
1727	chem_debug2, &
1728	chem_gasphase_on, &
1729	cs_heights, &
1730	cs_name, &
1731	cs_profile, &
1732	cs_profile_name, &
1733	cs_surface, &
1734	decycle_chem_lr, &
1735	decycle_chem_ns, &
1736	decycle_method, &
1737	emiss_factor_main, &
1738	emiss_factor_side, &
1739	icntrl, &
1740	main_street_id, &
1741	max_street_id, &
1742	my_steps, &
1743	nest_chemistry, &
1744	rcntrl, &
1745	side_street_id, &
1746	photolysis_scheme, &
1747	wall_csflux, &
1748	cs_vertical_gradient, &
1749	top_csflux, &
1750	surface_csflux, &
1751	surface_csflux_name, &
1752	cs_surface_initial_change, &
1753	cs_vertical_gradient_level, &
1754	! namelist parameters for emissions
1755	mode_emis, &
1756	time_fac_type, &
1757	daytype_mdh, &
1758	do_emis
1759
1760	!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
1761	!-- so this way we could prescribe a specific flux value for each species
1762	!> chemistry_parameters for initial profiles
1763	!> cs_names = 'O3', 'NO2', 'NO', ... to set initial profiles)
1764	!> cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
1765	!> cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
1766	!> cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, ..... (chem spcs conc at height lvls chem_heights(1,:)) etc.
1767	!> If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
1768	!> then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
1769
1770	!
1771	!-- Read chem namelist
1772
1773	INTEGER :: ier
1774	CHARACTER(LEN=64) :: text
1775	CHARACTER(LEN=8) :: solver_type
1776
1777	icntrl = 0
1778	rcntrl = 0.0_wp
1779	my_steps = 0.0_wp
1780	photolysis_scheme = 'simple'
1781	atol = 1.0_wp
1782	rtol = 0.01_wp
1783	!
1784	!-- Try to find chemistry package
1785	REWIND ( 11 )
1786	line = ' '
1787	DO WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
1788	READ ( 11, '(A)', END=20 ) line
1789	ENDDO
1790	BACKSPACE ( 11 )
1791	!
1792	!-- Read chemistry namelist
1793	READ ( 11, chemistry_parameters, ERR = 10, END = 20 )
1794	!
1795	!-- Enable chemistry model
1796	air_chemistry = .TRUE.
1797	GOTO 20
1798
1799	10 BACKSPACE( 11 )
1800	READ( 11 , '(A)') line
1801	CALL parin_fail_message( 'chemistry_parameters', line )
1802
1803	20 CONTINUE
1804
1805	!
1806	!-- check for emission mode for chem species
1807	IF ( (mode_emis /= 'PARAMETERIZED') .AND. ( mode_emis /= 'DEFAULT') .AND. (mode_emis /= 'PRE-PROCESSED' ) ) THEN
1808	message_string = 'Incorrect mode_emiss option select. Please check spelling'
1809	CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
1810	ENDIF
1811
1812	t_steps = my_steps
1813
1814	!-- Determine the number of user-defined profiles and append them to the
1815	!-- standard data output (data_output_pr)
1816	max_pr_cs_tmp = 0
1817	i = 1
1818
1819	DO WHILE ( data_output_pr(i) /= ' ' .AND. i <= 100 )
1820	IF ( TRIM( data_output_pr(i)(1:3)) == 'kc_' ) THEN
1821	max_pr_cs_tmp = max_pr_cs_tmp+1
1822	ENDIF
1823	i = i +1
1824	ENDDO
1825
1826	IF ( max_pr_cs_tmp > 0 ) THEN
1827	cs_pr_namelist_found = .TRUE.
1828	max_pr_cs = max_pr_cs_tmp
1829	ENDIF
1830
1831	! Set Solver Type
1832	IF(icntrl(3) == 0) THEN
1833	solver_type = 'rodas3' !Default
1834	ELSE IF(icntrl(3) == 1) THEN
1835	solver_type = 'ros2'
1836	ELSE IF(icntrl(3) == 2) THEN
1837	solver_type = 'ros3'
1838	ELSE IF(icntrl(3) == 3) THEN
1839	solver_type = 'ro4'
1840	ELSE IF(icntrl(3) == 4) THEN
1841	solver_type = 'rodas3'
1842	ELSE IF(icntrl(3) == 5) THEN
1843	solver_type = 'rodas4'
1844	ELSE IF(icntrl(3) == 6) THEN
1845	solver_type = 'Rang3'
1846	ELSE
1847	IF(myid == 0) write(0,*) 'illegal solver type ' !kk chane to PALM error message
1848	call MPI_Abort (MPI_COMM_WORLD, 1, ier)
1849	END IF
1850
1851	write(text,*) 'gas_phase chemistry: solver_type = ',trim(solver_type)
1852	!kk Has to be changed to right calling sequence
1853	!kk CALL location_message( trim(text), .FALSE. )
1854	IF(myid == 0) THEN
1855	write(9,*) ' '
1856	write(9,*) 'kpp setup '
1857	write(9,*) ' '
1858	write(9,*) ' gas_phase chemistry: solver_type = ',trim(solver_type)
1859	write(9,*) ' '
1860	write(9,*) ' Hstart = ',rcntrl(3)
1861	write(9,*) ' FacMin = ',rcntrl(4)
1862	write(9,*) ' FacMax = ',rcntrl(5)
1863	write(9,*) ' '
1864	IF(vl_dim > 1) THEN
1865	write(9,*) ' Vector mode vektor length = ',vl_dim
1866	ELSE
1867	write(9,*) ' Scalar mode'
1868	ENDIF
1869	write(9,*) ' '
1870	END IF
1871
1872	RETURN
1873
1874	END SUBROUTINE chem_parin
1875
1876	!
1877	!------------------------------------------------------------------------------!
1878	!
1879	! Description:
1880	! ------------
1881	!> Subroutine calculating prognostic equations for chemical species
1882	!> (vector-optimized).
1883	!> Routine is called separately for each chemical species over a loop from
1884	!> prognostic_equations.
1885	!------------------------------------------------------------------------------!
1886	SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m, &
1887	pr_init_cs, ilsp )
1888
1889	USE advec_s_pw_mod, &
1890	ONLY: advec_s_pw
1891	USE advec_s_up_mod, &
1892	ONLY: advec_s_up
1893	USE advec_ws, &
1894	ONLY: advec_s_ws
1895	USE diffusion_s_mod, &
1896	ONLY: diffusion_s
1897	USE indices, &
1898	ONLY: nxl, nxr, nyn, nys, wall_flags_0
1899	USE pegrid
1900	USE surface_mod, &
1901	ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
1902	surf_usm_v
1903
1904	IMPLICIT NONE
1905
1906	INTEGER :: i !< running index
1907	INTEGER :: j !< running index
1908	INTEGER :: k !< running index
1909
1910	INTEGER(iwp),INTENT(IN) :: ilsp !<
1911
1912	REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !<
1913
1914	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar !<
1915	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p !<
1916	REAL(wp), DIMENSION(:,:,:), POINTER :: tcs_scalar_m !<
1917
1918
1919	!
1920	!-- Tendency terms for chemical species
1921	tend = 0.0_wp
1922	!
1923	!-- Advection terms
1924	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1925	IF ( ws_scheme_sca ) THEN
1926	CALL advec_s_ws( cs_scalar, 'kc' )
1927	ELSE
1928	CALL advec_s_pw( cs_scalar )
1929	ENDIF
1930	ELSE
1931	CALL advec_s_up( cs_scalar )
1932	ENDIF
1933	!
1934	!-- Diffusion terms (the last three arguments are zero)
1935	CALL diffusion_s( cs_scalar, &
1936	surf_def_h(0)%cssws(ilsp,:), &
1937	surf_def_h(1)%cssws(ilsp,:), &
1938	surf_def_h(2)%cssws(ilsp,:), &
1939	surf_lsm_h%cssws(ilsp,:), &
1940	surf_usm_h%cssws(ilsp,:), &
1941	surf_def_v(0)%cssws(ilsp,:), &
1942	surf_def_v(1)%cssws(ilsp,:), &
1943	surf_def_v(2)%cssws(ilsp,:), &
1944	surf_def_v(3)%cssws(ilsp,:), &
1945	surf_lsm_v(0)%cssws(ilsp,:), &
1946	surf_lsm_v(1)%cssws(ilsp,:), &
1947	surf_lsm_v(2)%cssws(ilsp,:), &
1948	surf_lsm_v(3)%cssws(ilsp,:), &
1949	surf_usm_v(0)%cssws(ilsp,:), &
1950	surf_usm_v(1)%cssws(ilsp,:), &
1951	surf_usm_v(2)%cssws(ilsp,:), &
1952	surf_usm_v(3)%cssws(ilsp,:) )
1953	!
1954	!-- Prognostic equation for chemical species
1955	DO i = nxl, nxr
1956	DO j = nys, nyn
1957	DO k = nzb+1, nzt
1958	cs_scalar_p(k,j,i) = cs_scalar(k,j,i) &
1959	+ ( dt_3d * &
1960	( tsc(2) * tend(k,j,i) &
1961	+ tsc(3) * tcs_scalar_m(k,j,i) &
1962	) &
1963	- tsc(5) * rdf_sc(k) &
1964	* ( cs_scalar(k,j,i) - pr_init_cs(k) ) &
1965	) &
1966	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1967
1968	IF ( cs_scalar_p(k,j,i) < 0.0_wp ) cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)
1969	ENDDO
1970	ENDDO
1971	ENDDO
1972	!
1973	!-- Calculate tendencies for the next Runge-Kutta step
1974	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1975	IF ( intermediate_timestep_count == 1 ) THEN
1976	DO i = nxl, nxr
1977	DO j = nys, nyn
1978	DO k = nzb+1, nzt
1979	tcs_scalar_m(k,j,i) = tend(k,j,i)
1980	ENDDO
1981	ENDDO
1982	ENDDO
1983	ELSEIF ( intermediate_timestep_count < &
1984	intermediate_timestep_count_max ) THEN
1985	DO i = nxl, nxr
1986	DO j = nys, nyn
1987	DO k = nzb+1, nzt
1988	tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i) &
1989	+ 5.3125_wp * tcs_scalar_m(k,j,i)
1990	ENDDO
1991	ENDDO
1992	ENDDO
1993	ENDIF
1994	ENDIF
1995
1996	END SUBROUTINE chem_prognostic_equations
1997
1998	!------------------------------------------------------------------------------!
1999	!
2000	! Description:
2001	! ------------
2002	!> Subroutine calculating prognostic equations for chemical species
2003	!> (cache-optimized).
2004	!> Routine is called separately for each chemical species over a loop from
2005	!> prognostic_equations.
2006	!------------------------------------------------------------------------------!
2007	SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs, &
2008	i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, &
2009	flux_l_cs, diss_l_cs )
2010	USE pegrid
2011	USE advec_ws, ONLY: advec_s_ws
2012	USE advec_s_pw_mod, ONLY: advec_s_pw
2013	USE advec_s_up_mod, ONLY: advec_s_up
2014	USE diffusion_s_mod, ONLY: diffusion_s
2015	USE indices, ONLY: wall_flags_0
2016	USE surface_mod, ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
2017	surf_usm_v
2018
2019
2020	IMPLICIT NONE
2021
2022	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p, cs_scalar, tcs_scalar_m
2023
2024	INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp
2025	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: flux_s_cs !<
2026	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: diss_s_cs !<
2027	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: flux_l_cs !<
2028	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: diss_l_cs !<
2029	REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !<
2030
2031	!-- local variables
2032
2033	INTEGER :: k
2034	!
2035	!-- Tendency-terms for chem spcs.
2036	tend(:,j,i) = 0.0_wp
2037	!
2038	!-- Advection terms
2039	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2040	IF ( ws_scheme_sca ) THEN
2041	CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs, &
2042	flux_l_cs, diss_l_cs, i_omp_start, tn )
2043	ELSE
2044	CALL advec_s_pw( i, j, cs_scalar )
2045	ENDIF
2046	ELSE
2047	CALL advec_s_up( i, j, cs_scalar )
2048	ENDIF
2049
2050	!
2051
2052	!-- Diffusion terms (the last three arguments are zero)
2053
2054	CALL diffusion_s( i, j, cs_scalar, &
2055	surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:), &
2056	surf_def_h(2)%cssws(ilsp,:), &
2057	surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:), &
2058	surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:), &
2059	surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:), &
2060	surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:), &
2061	surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:), &
2062	surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:), &
2063	surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
2064
2065	!
2066	!-- Prognostic equation for chem spcs
2067	DO k = nzb+1, nzt
2068	cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d * &
2069	( tsc(2) * tend(k,j,i) + &
2070	tsc(3) * tcs_scalar_m(k,j,i) ) &
2071	- tsc(5) * rdf_sc(k) &
2072	* ( cs_scalar(k,j,i) - pr_init_cs(k) ) &
2073	) &
2074	* MERGE( 1.0_wp, 0.0_wp, &
2075	BTEST( wall_flags_0(k,j,i), 0 ) &
2076	)
2077
2078	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
2079	ENDDO
2080
2081	!
2082	!-- Calculate tendencies for the next Runge-Kutta step
2083	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2084	IF ( intermediate_timestep_count == 1 ) THEN
2085	DO k = nzb+1, nzt
2086	tcs_scalar_m(k,j,i) = tend(k,j,i)
2087	ENDDO
2088	ELSEIF ( intermediate_timestep_count < &
2089	intermediate_timestep_count_max ) THEN
2090	DO k = nzb+1, nzt
2091	tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
2092	5.3125_wp * tcs_scalar_m(k,j,i)
2093	ENDDO
2094	ENDIF
2095	ENDIF
2096
2097	END SUBROUTINE chem_prognostic_equations_ij
2098
2099	!
2100	!------------------------------------------------------------------------------!
2101	!
2102	! Description:
2103	! ------------
2104	!> Subroutine to read restart data of chemical species
2105	!------------------------------------------------------------------------------!
2106
2107	SUBROUTINE chem_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
2108	nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc, &
2109	nys_on_file, tmp_3d, found )
2110
2111	USE control_parameters
2112
2113	USE indices
2114
2115	USE pegrid
2116
2117	IMPLICIT NONE
2118
2119	CHARACTER (LEN=20) :: spc_name_av !<
2120
2121	INTEGER(iwp) :: i, lsp !<
2122	INTEGER(iwp) :: k !<
2123	INTEGER(iwp) :: nxlc !<
2124	INTEGER(iwp) :: nxlf !<
2125	INTEGER(iwp) :: nxl_on_file !<
2126	INTEGER(iwp) :: nxrc !<
2127	INTEGER(iwp) :: nxrf !<
2128	INTEGER(iwp) :: nxr_on_file !<
2129	INTEGER(iwp) :: nync !<
2130	INTEGER(iwp) :: nynf !<
2131	INTEGER(iwp) :: nyn_on_file !<
2132	INTEGER(iwp) :: nysc !<
2133	INTEGER(iwp) :: nysf !<
2134	INTEGER(iwp) :: nys_on_file !<
2135
2136	LOGICAL, INTENT(OUT) :: found
2137
2138	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
2139
2140
2141	found = .FALSE.
2142
2143
2144	IF ( ALLOCATED(chem_species) ) THEN
2145
2146	DO lsp = 1, nspec
2147
2148	!< for time-averaged chemical conc.
2149	spc_name_av = TRIM(chem_species(lsp)%name)//'_av'
2150
2151	IF (restart_string(1:length) == TRIM(chem_species(lsp)%name) ) &
2152	THEN
2153	!< read data into tmp_3d
2154	IF ( k == 1 ) READ ( 13 ) tmp_3d
2155	!< fill ..%conc in the restart run
2156	chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp, &
2157	nxlc-nbgp:nxrc+nbgp) = &
2158	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2159	found = .TRUE.
2160	ELSEIF (restart_string(1:length) == spc_name_av ) THEN
2161	IF ( k == 1 ) READ ( 13 ) tmp_3d
2162	chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp, &
2163	nxlc-nbgp:nxrc+nbgp) = &
2164	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2165	found = .TRUE.
2166	ENDIF
2167
2168	ENDDO
2169
2170	ENDIF
2171
2172
2173	END SUBROUTINE chem_rrd_local
2174	!
2175	!-------------------------------------------------------------------------------!
2176	!> Description:
2177	!> Calculation of horizontally averaged profiles
2178	!> This routine is called for every statistic region (sr) defined by the user,
2179	!> but at least for the region "total domain" (sr=0).
2180	!> quantities.
2181	!------------------------------------------------------------------------------!
2182	SUBROUTINE chem_statistics( mode, sr, tn )
2183
2184
2185	USE arrays_3d
2186	USE indices
2187	USE kinds
2188	USE pegrid
2189	USE statistics
2190
2191	USE user
2192
2193	IMPLICIT NONE
2194
2195	CHARACTER (LEN=*) :: mode !<
2196
2197	INTEGER(iwp) :: i !< running index on x-axis
2198	INTEGER(iwp) :: j !< running index on y-axis
2199	INTEGER(iwp) :: k !< vertical index counter
2200	INTEGER(iwp) :: sr !< statistical region
2201	INTEGER(iwp) :: tn !< thread number
2202	INTEGER(iwp) :: n !<
2203	INTEGER(iwp) :: m !<
2204	INTEGER(iwp) :: lpr !< running index chem spcs
2205	! REAL(wp), &
2206	! DIMENSION(dots_num_palm+1:dots_max) :: &
2207	! ts_value_l !<
2208
2209	IF ( mode == 'profiles' ) THEN
2210	!
2211	!-- Sample on how to calculate horizontally averaged profiles of user-
2212	!-- defined quantities. Each quantity is identified by the index
2213	!-- "pr_palm+#" where "#" is an integer starting from 1. These
2214	!-- user-profile-numbers must also be assigned to the respective strings
2215	!-- given by data_output_pr_cs in routine user_check_data_output_pr.
2216	!-- hom(:,:,:,:) = dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
2217	! wpt), dim-4 = statistical region.
2218
2219	!$OMP DO
2220	DO i = nxl, nxr
2221	DO j = nys, nyn
2222	DO k = nzb, nzt+1
2223	DO lpr = 1, cs_pr_count
2224
2225	sums_l(k,pr_palm+lpr,tn) = sums_l(k,pr_palm+lpr,tn) + &
2226	chem_species(cs_pr_index(lpr))%conc(k,j,i) * &
2227	rmask(j,i,sr) * &
2228	MERGE( 1.0_wp, 0.0_wp, &
2229	BTEST( wall_flags_0(k,j,i), 22 ) )
2230	ENDDO
2231	ENDDO
2232	ENDDO
2233	ENDDO
2234	ELSEIF ( mode == 'time_series' ) THEN
2235	CALL location_message( 'Time series not calculated for chemistry', .TRUE. )
2236	ENDIF
2237
2238	END SUBROUTINE chem_statistics
2239	!
2240	!------------------------------------------------------------------------------!
2241	!
2242	! Description:
2243	! ------------
2244	!> Subroutine for swapping of timelevels for chemical species
2245	!> called out from subroutine swap_timelevel
2246	!------------------------------------------------------------------------------!
2247
2248	SUBROUTINE chem_swap_timelevel (level)
2249
2250	IMPLICIT none
2251
2252	INTEGER,INTENT(IN) :: level
2253	!-- local variables
2254	INTEGER :: lsp
2255
2256
2257	IF ( level == 0 ) THEN
2258	DO lsp=1, nvar
2259	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp)
2260	chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp)
2261	ENDDO
2262	ELSE
2263	DO lsp=1, nvar
2264	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp)
2265	chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp)
2266	ENDDO
2267	ENDIF
2268
2269	RETURN
2270	END SUBROUTINE chem_swap_timelevel
2271	!
2272	!------------------------------------------------------------------------------!
2273	!
2274	! Description:
2275	! ------------
2276	!> Subroutine to write restart data for chemistry model
2277	!------------------------------------------------------------------------------!
2278	SUBROUTINE chem_wrd_local
2279
2280	IMPLICIT NONE
2281
2282	INTEGER(iwp) :: lsp !<
2283
2284	DO lsp = 1, nspec
2285	CALL wrd_write_string( TRIM( chem_species(lsp)%name ))
2286	WRITE ( 14 ) chem_species(lsp)%conc
2287	CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
2288	WRITE ( 14 ) chem_species(lsp)%conc_av
2289	ENDDO
2290
2291	END SUBROUTINE chem_wrd_local
2292	END MODULE chemistry_model_mod
2293

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |