Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |