Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |