Home

Context Navigation

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

Last change on this file since 3435 was 3435, checked in by gronemeier, 6 years ago
new: terrain-following masked output; bugfixes: increase vertical dimension of gamma_w_green_sat by 1, add checks for masked output for chemistry_model and radiation_model, reordered calls to xxx_define_netcdf_grid in masked output part
Property svn:keywords set to `Id`
File size: 98.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |