Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |