Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |