Home

Context Navigation

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

Last change on this file since 2766 was 2766, checked in by kanani, 6 years ago
Removal of chem directive, plus minor changes
Property svn:keywords set to `Id`
File size: 72.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |