Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |