Home

Context Navigation

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

Last change on this file since 3025 was 3014, checked in by maronga, 7 years ago
series of bugfixes
Property svn:keywords set to `Id`
File size: 75.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |