Home

Context Navigation

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

Last change on this file since 3246 was 3246, checked in by sward, 6 years ago
Added error handling for wrong input parameters
Property svn:keywords set to `Id`
File size: 75.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |