Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |