Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |