Home

Context Navigation

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

Last change on this file since 3058 was 3045, checked in by Giersch, 7 years ago
Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled
Property svn:keywords set to `Id`
File size: 75.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |