Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |