Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |