Home

Context Navigation

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

Last change on this file since 3376 was 3373, checked in by kanani, 6 years ago
Fix cpp directives and error messages
Property svn:keywords set to `Id`
File size: 95.5 KB

Rev	Line
[2425]	1	!> @file chemistry_model_mod.f90
	2	!------------------------------------------------------------------------------!
[2828]	3	! This file is part of the PALM model system.
[2425]	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	!
[3298]	17	! Copyright 2017-2018 Leibniz Universitaet Hannover
[2828]	18	! Copyright 2017-2018 Karlsruhe Institute of Technology
	19	! Copyright 2017-2018 Freie Universitaet Berlin
[2425]	20	!------------------------------------------------------------------------------!
	21	!
	22	! Current revisions:
	23	! -----------------
[3298]	24	!
	25	!
	26	! Former revisions:
	27	! -----------------
	28	! $Id: chemistry_model_mod.f90 3373 2018-10-18 15:25:56Z suehring $
[3373]	29	! Remove MPI_Abort, replace by message
	30	!
	31	! 3318 2018-10-08 11:43:01Z sward
[3318]	32	! Fixed faulty syntax of message string
	33	!
	34	! 3298 2018-10-02 12:21:11Z kanani
[3298]	35	! Add remarks (kanani)
[3292]	36	! Merge with trunk, replaced cloud_physics by bulk_cloud_model 28.09.2018 forkel
[3281]	37	! Subroutines header and chem_check_parameters added 25.09.2018 basit
[3190]	38	! Removed chem_emission routine now declared in chem_emissions.f90 30.07.2018 ERUSSO
	39	! Introduced emissions namelist parameters 30.07.2018 ERUSSO
[3188]	40	!
	41	! Timestep steering added in subroutine chem_integrate_ij and
	42	! output of chosen solver in chem_parin added 30.07.2018 ketelsen
	43	!
	44	! chem_check_data_output_pr: added unit for PM compounds 20.07.2018 forkel
	45	! replaced : by nzb+1:nzt for pt,q,ql (found by kk) 18.07.2018 forkel
[3140]	46	! debugged restart run for chem species 06.07.2018 basit
	47	! reorganized subroutines in alphabetical order. 27.06.2018 basit
	48	! subroutine chem_parin updated for profile output 27.06.2018 basit
[3114]	49	! Added humidity arrays to USE section and tmp_qvap in chem_integrate 26.6.2018 forkel
	50	! Merged chemistry with with trunk (nzb_do and nzt_do in 3d output) 26.6.2018 forkel
	51	!
[3140]	52	! reorganized subroutines in alphabetical order. basit 22.06.2018
	53	! subroutine chem_parin updated for profile output basit 22.06.2018
	54	! subroutine chem_statistics added
	55	! subroutine chem_check_data_output_pr add 21.06.2018 basit
	56	! subroutine chem_data_output_mask added 20.05.2018 basit
	57	! subroutine chem_data_output_2d added 20.05.2018 basit
	58	! subroutine chem_statistics added 04.06.2018 basit
[2914]	59	! subroutine chem_emissions: Set cssws to zero before setting values 20.03.2018 forkel
[2900]	60	! subroutine chem_emissions: Introduced different conversion factors
	61	! for PM and gaseous compounds 15.03.2018 forkel
[2888]	62	! subroutine chem_emissions updated to take variable number of chem_spcs and
[3140]	63	! emission factors. 13.03.2018 basit
	64	! chem_boundary_conds_decycle improved. 05.03.2018 basit
	65	! chem_boundary_conds_decycle subroutine added 21.02.2018 basit
[3188]	66	! chem_init_profiles subroutines re-activated after correction 21.02.2018 basit
[2616]	67	!
[2828]	68	!
[3298]	69	! 3293 2018-09-28 12:45:20Z forkel
[3287]	70	! Modularization of all bulk cloud physics code components
	71	!
	72	! 3248 2018-09-14 09:42:06Z sward
	73	! Minor formating changes
	74	!
	75	! 3246 2018-09-13 15:14:50Z sward
	76	! Added error handling for input namelist via parin_fail_message
	77	!
	78	! 3241 2018-09-12 15:02:00Z raasch
[3228]	79	! +nest_chemistry
	80	!
	81	! 3209 2018-08-27 16:58:37Z suehring
	82	! Rename flags indicating outflow boundary conditions
	83	!
	84	! 3182 2018-07-27 13:36:03Z suehring
	85	! Revise output of surface quantities in case of overhanging structures
	86	!
	87	! 3045 2018-05-28 07:55:41Z Giersch
[3114]	88	! error messages revised
	89	!
	90	! 3014 2018-05-09 08:42:38Z maronga
	91	! Bugfix: nzb_do and nzt_do were not used for 3d data output
	92	!
	93	! 3004 2018-04-27 12:33:25Z Giersch
	94	! Comment concerning averaged data output added
	95	!
	96	! 2932 2018-03-26 09:39:22Z maronga
[2981]	97	! renamed chemistry_par to chemistry_parameters
	98	!
	99	! 2894 2018-03-15 09:17:58Z Giersch
[2914]	100	! Calculations of the index range of the subdomain on file which overlaps with
	101	! the current subdomain are already done in read_restart_data_mod,
	102	! chem_last_actions was renamed to chem_wrd_local, chem_read_restart_data was
	103	! renamed to chem_rrd_local, chem_write_var_list was renamed to
	104	! chem_wrd_global, chem_read_var_list was renamed to chem_rrd_global,
	105	! chem_skip_var_list has been removed, variable named found has been
	106	! introduced for checking if restart data was found, reading of restart strings
	107	! has been moved completely to read_restart_data_mod, chem_rrd_local is already
	108	! inside the overlap loop programmed in read_restart_data_mod, todo list has
	109	! bees extended, redundant characters in chem_wrd_local have been removed,
	110	! the marker * end chemistry * is not necessary anymore, strings and their
	111	! respective lengths are written out and read now in case of restart runs to
	112	! get rid of prescribed character lengths
	113	!
	114	! 2815 2018-02-19 11:29:57Z suehring
[2828]	115	! Bugfix in restart mechanism,
	116	! rename chem_tendency to chem_prognostic_equations,
	117	! implement vector-optimized version of chem_prognostic_equations,
	118	! some clean up (incl. todo list)
[2656]	119	!
[2828]	120	! 2773 2018-01-30 14:12:54Z suehring
	121	! Declare variables required for nesting as public
[2656]	122	!
[2828]	123	! 2772 2018-01-29 13:10:35Z suehring
	124	! Bugfix in string handling
[2635]	125	!
[2828]	126	! 2768 2018-01-24 15:38:29Z kanani
	127	! Shorten lines to maximum length of 132 characters
[2633]	128	!
[2828]	129	! 2766 2018-01-22 17:17:47Z kanani
	130	! Removed preprocessor directive __chem
[2633]	131	!
[2828]	132	! 2756 2018-01-16 18:11:14Z suehring
	133	! Fill values in 3D output introduced.
[2627]	134	!
[2828]	135	! 2718 2018-01-02 08:49:38Z maronga
	136	! Initial revision
	137	!
[2452]	138	!
[2828]	139	!
	140	!
	141	! Authors:
	142	! --------
	143	! @author Renate Forkel
	144	! @author Farah Kanani-Suehring
	145	! @author Klaus Ketelsen
	146	! @author Basit Khan
	147	!
	148	!
[2426]	149	!------------------------------------------------------------------------------!
[2425]	150	! Description:
	151	! ------------
[2592]	152	!> Chemistry model for PALM-4U
[2914]	153	!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
	154	!> allowed to use the chemistry model in a precursor run and additionally
	155	!> not using it in a main run
[2828]	156	!> @todo Update/clean-up todo list! (FK)
	157	!> @todo Set proper fill values (/= 0) for chem output arrays! (FK)
	158	!> @todo Add routine chem_check_parameters, add checks for inconsistent or
	159	!> unallowed parameter settings.
	160	!> CALL of chem_check_parameters from check_parameters. (FK)
	161	!> @todo Make routine chem_header available, CALL from header.f90
	162	!> (see e.g. how it is done in routine lsm_header in
	163	!> land_surface_model_mod.f90). chem_header should include all setup
	164	!> info about chemistry parameter settings. (FK)
[2592]	165	!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
	166	!> @todo Separate boundary conditions for each chem spcs to be implemented
	167	!> @todo Currently only total concentration are calculated. Resolved, parameterized
	168	!> and chemistry fluxes although partially and some completely coded but
	169	!> are not operational/activated in this version. bK.
	170	!> @todo slight differences in passive scalar and chem spcs when chem reactions
	171	!> turned off. Need to be fixed. bK
[2828]	172	!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
[2592]	173	!> @todo chemistry error messages
[2425]	174	!> @todo Format this module according to PALM coding standard (see e.g. module
	175	!> template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or
	176	!> D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16)
[2482]	177	!
[2425]	178	!------------------------------------------------------------------------------!
	179
[3281]	180	MODULE chemistry_model_mod
[2828]	181
[3281]	182	USE kinds, ONLY: wp, iwp
[3282]	183	USE indices, ONLY: nz, nzb,nzt,nysg,nyng,nxlg,nxrg, &
[3114]	184	nys,nyn,nx,nxl,nxr,ny,wall_flags_0
[3281]	185	USE pegrid, ONLY: myid, threads_per_task
[3293]	186
	187	USE bulk_cloud_model_mod, &
	188	ONLY: bulk_cloud_model
	189
	190	USE control_parameters, ONLY: bc_lr_cyc, bc_ns_cyc, dt_3d, humidity, &
[3282]	191	initializing_actions, message_string, &
	192	omega, tsc, intermediate_timestep_count, &
	193	intermediate_timestep_count_max, &
[3287]	194	timestep_scheme, use_prescribed_profile_data, ws_scheme_sca
	195	USE arrays_3d, ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu
[3282]	196	USE chem_gasphase_mod, ONLY: nspec, spc_names, nkppctrl, nmaxfixsteps, &
	197	t_steps, chem_gasphase_integrate, vl_dim, &
	198	nvar, nreact, atol, rtol, nphot, phot_names
[3281]	199	USE cpulog, ONLY: cpu_log, log_point
[2425]	200
[3281]	201	USE chem_modules
[2615]	202
[3281]	203	USE statistics
[2615]	204
[3281]	205	IMPLICIT NONE
	206	PRIVATE
	207	SAVE
[2425]	208
[3281]	209	LOGICAL :: nest_chemistry = .TRUE. !< flag for nesting mode of chemical species, independent on parent or not
[3228]	210	!
[2425]	211	!- Define chemical variables
[3281]	212	TYPE species_def
	213	CHARACTER(LEN=8) :: name
	214	CHARACTER(LEN=16) :: unit
	215	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc
	216	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_av
	217	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: conc_p
	218	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: tconc_m
	219	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: cssws_av
	220	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: flux_s_cs
	221	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:) :: diss_s_cs
	222	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: flux_l_cs
	223	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:) :: diss_l_cs
	224	REAL(kind=wp),ALLOCATABLE,DIMENSION(:) :: conc_pr_init
	225	END TYPE species_def
[2425]	226
[2592]	227
[3281]	228	TYPE photols_def
	229	CHARACTER(LEN=8) :: name
	230	CHARACTER(LEN=16) :: unit
	231	REAL(kind=wp),POINTER,DIMENSION(:,:,:) :: freq
	232	END TYPE photols_def
[2425]	233
	234
[3281]	235	PUBLIC species_def
	236	PUBLIC photols_def
[2425]	237
[3228]	238
[3281]	239	TYPE(species_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: chem_species
	240	TYPE(photols_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC :: phot_frequen
[2425]	241
[3281]	242	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_1
	243	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_2
	244	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_3
	245	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: spec_conc_av
	246	REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET :: freq_1
[2425]	247
[3281]	248	INTEGER,DIMENSION(nkppctrl) :: icntrl ! Fine tuning kpp
	249	REAL(kind=wp),DIMENSION(nkppctrl) :: rcntrl ! Fine tuning kpp
	250	LOGICAL :: decycle_chem_lr = .FALSE.
	251	LOGICAL :: decycle_chem_ns = .FALSE.
	252	CHARACTER (LEN=20), DIMENSION(4) :: decycle_method = &
[2854]	253	(/'dirichlet','dirichlet','dirichlet','dirichlet'/)
	254	!< Decycling method at horizontal boundaries,
	255	!< 1=left, 2=right, 3=south, 4=north
	256	!< dirichlet = initial size distribution and
	257	!< chemical composition set for the ghost and
	258	!< first three layers
	259	!< neumann = zero gradient
[2425]	260
[3282]	261	REAL(kind=wp), PUBLIC :: cs_time_step
[3281]	262	CHARACTER(10), PUBLIC :: photolysis_scheme
[3114]	263	! 'constant',
[2592]	264	! 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180
	265	! 'fastj' (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
[2535]	266
[3281]	267	PUBLIC nest_chemistry
	268	PUBLIC nspec
[3282]	269	PUBLIC nreact
[3281]	270	PUBLIC nvar
	271	PUBLIC spc_names
	272	PUBLIC spec_conc_2
[2425]	273
[3281]	274	!- Interface section
	275	INTERFACE chem_3d_data_averaging
	276	MODULE PROCEDURE chem_3d_data_averaging
	277	END INTERFACE chem_3d_data_averaging
[3228]	278
[3281]	279	INTERFACE chem_boundary_conds
	280	MODULE PROCEDURE chem_boundary_conds
	281	MODULE PROCEDURE chem_boundary_conds_decycle
	282	END INTERFACE chem_boundary_conds
[2425]	283
[3281]	284	INTERFACE chem_check_data_output
	285	MODULE PROCEDURE chem_check_data_output
	286	END INTERFACE chem_check_data_output
[2425]	287
[3281]	288	INTERFACE chem_data_output_2d
	289	MODULE PROCEDURE chem_data_output_2d
	290	END INTERFACE chem_data_output_2d
[2425]	291
[3281]	292	INTERFACE chem_data_output_3d
	293	MODULE PROCEDURE chem_data_output_3d
	294	END INTERFACE chem_data_output_3d
[2425]	295
[3281]	296	INTERFACE chem_data_output_mask
	297	MODULE PROCEDURE chem_data_output_mask
	298	END INTERFACE chem_data_output_mask
[3085]	299
[3281]	300	INTERFACE chem_check_data_output_pr
	301	MODULE PROCEDURE chem_check_data_output_pr
	302	END INTERFACE chem_check_data_output_pr
[3085]	303
[3281]	304	INTERFACE chem_check_parameters
	305	MODULE PROCEDURE chem_check_parameters
	306	END INTERFACE chem_check_parameters
[2425]	307
[3281]	308	INTERFACE chem_define_netcdf_grid
	309	MODULE PROCEDURE chem_define_netcdf_grid
	310	END INTERFACE chem_define_netcdf_grid
[2425]	311
[3281]	312	INTERFACE chem_header
	313	MODULE PROCEDURE chem_header
	314	END INTERFACE chem_header
[2425]	315
[3281]	316	INTERFACE chem_init
	317	MODULE PROCEDURE chem_init
	318	END INTERFACE chem_init
[2425]	319
[3281]	320	INTERFACE chem_init_profiles
	321	MODULE PROCEDURE chem_init_profiles
	322	END INTERFACE chem_init_profiles
[2425]	323
[3281]	324	INTERFACE chem_integrate
	325	MODULE PROCEDURE chem_integrate_ij
	326	END INTERFACE chem_integrate
[2425]	327
[3281]	328	INTERFACE chem_parin
	329	MODULE PROCEDURE chem_parin
	330	END INTERFACE chem_parin
[2425]	331
[3281]	332	INTERFACE chem_prognostic_equations
	333	MODULE PROCEDURE chem_prognostic_equations
	334	MODULE PROCEDURE chem_prognostic_equations_ij
	335	END INTERFACE chem_prognostic_equations
[3228]	336
[3281]	337	INTERFACE chem_rrd_local
	338	MODULE PROCEDURE chem_rrd_local
	339	END INTERFACE chem_rrd_local
[2467]	340
[3281]	341	INTERFACE chem_statistics
	342	MODULE PROCEDURE chem_statistics
	343	END INTERFACE chem_statistics
[3085]	344
[3281]	345	INTERFACE chem_swap_timelevel
	346	MODULE PROCEDURE chem_swap_timelevel
	347	END INTERFACE chem_swap_timelevel
	348
	349	INTERFACE chem_wrd_local
	350	MODULE PROCEDURE chem_wrd_local
	351	END INTERFACE chem_wrd_local
[2482]	352
[2615]	353
[3281]	354	PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
	355	chem_check_data_output_pr, chem_check_parameters, &
	356	chem_data_output_2d, chem_data_output_3d, chem_data_output_mask, &
	357	chem_define_netcdf_grid, chem_header, chem_init, &
	358	chem_init_profiles, chem_integrate, chem_parin, &
	359	chem_prognostic_equations, chem_rrd_local, &
	360	chem_statistics, chem_swap_timelevel, chem_wrd_local
[2425]	361
[3281]	362
	363
[2425]	364	CONTAINS
	365
[3228]	366	!
[2425]	367	!------------------------------------------------------------------------------!
[2535]	368	!
[2425]	369	! Description:
	370	! ------------
[3228]	371	!> Subroutine for averaging 3D data of chemical species. Due to the fact that
	372	!> the averaged chem arrays are allocated in chem_init, no if-query concerning
	373	!> the allocation is required (in any mode). Attention: If you just specify an
	374	!> averaged output quantity in the _p3dr file during restarts the first output
	375	!> includes the time between the beginning of the restart run and the first
	376	!> output time (not necessarily the whole averaging_interval you have
	377	!> specified in your _p3d/_p3dr file )
	378	!------------------------------------------------------------------------------!
	379
[3281]	380	SUBROUTINE chem_3d_data_averaging ( mode, variable )
[3228]	381
[3281]	382	USE control_parameters
[3228]	383
[3281]	384	USE indices
[3228]	385
[3281]	386	USE kinds
[3228]	387
[3281]	388	USE surface_mod, &
	389	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
[3228]	390
[3281]	391	IMPLICIT NONE
[3228]	392
[3281]	393	CHARACTER (LEN=*) :: mode !<
	394	CHARACTER (LEN=*) :: variable !<
[3228]	395
[3281]	396	LOGICAL :: match_def !< flag indicating natural-type surface
	397	LOGICAL :: match_lsm !< flag indicating natural-type surface
	398	LOGICAL :: match_usm !< flag indicating urban-type surface
[3228]	399
[3281]	400	INTEGER(iwp) :: i !< grid index x direction
	401	INTEGER(iwp) :: j !< grid index y direction
	402	INTEGER(iwp) :: k !< grid index z direction
	403	INTEGER(iwp) :: m !< running index surface type
	404	INTEGER(iwp) :: lsp !< running index for chem spcs
[3228]	405
[3281]	406
	407	IF ( mode == 'allocate' ) THEN
	408	DO lsp = 1, nspec
	409	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
	410	chem_species(lsp)%conc_av = 0.0_wp
	411
	412	ENDIF
	413	ENDDO
	414
	415	ELSEIF ( mode == 'sum' ) THEN
	416
	417
	418	DO lsp = 1, nspec
	419	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
	420	DO i = nxlg, nxrg
	421	DO j = nysg, nyng
	422	DO k = nzb, nzt+1
	423	chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + &
	424	chem_species(lsp)%conc(k,j,i)
	425	ENDDO
[3228]	426	ENDDO
	427	ENDDO
[3281]	428	ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') ) THEN
	429	DO i = nxl, nxr
	430	DO j = nys, nyn
	431	match_def = surf_def_h(0)%start_index(j,i) <= &
	432	surf_def_h(0)%end_index(j,i)
	433	match_lsm = surf_lsm_h%start_index(j,i) <= &
	434	surf_lsm_h%end_index(j,i)
	435	match_usm = surf_usm_h%start_index(j,i) <= &
	436	surf_usm_h%end_index(j,i)
[3228]	437
[3281]	438	IF ( match_def ) THEN
	439	m = surf_def_h(0)%end_index(j,i)
	440	chem_species(lsp)%cssws_av(j,i) = &
	441	chem_species(lsp)%cssws_av(j,i) + &
	442	surf_def_h(0)%cssws(lsp,m)
	443	ELSEIF ( match_lsm .AND. .NOT. match_usm ) THEN
	444	m = surf_lsm_h%end_index(j,i)
	445	chem_species(lsp)%cssws_av(j,i) = &
	446	chem_species(lsp)%cssws_av(j,i) + &
	447	surf_lsm_h%cssws(lsp,m)
	448	ELSEIF ( match_usm ) THEN
	449	m = surf_usm_h%end_index(j,i)
	450	chem_species(lsp)%cssws_av(j,i) = &
	451	chem_species(lsp)%cssws_av(j,i) + &
	452	surf_usm_h%cssws(lsp,m)
	453	ENDIF
	454	ENDDO
[3228]	455	ENDDO
[3281]	456	ENDIF
	457	ENDDO
[3228]	458
[3281]	459	ELSEIF ( mode == 'average' ) THEN
	460	DO lsp = 1, nspec
	461	IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
	462	DO i = nxlg, nxrg
	463	DO j = nysg, nyng
	464	DO k = nzb, nzt+1
	465	chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	466	ENDDO
[3228]	467	ENDDO
	468	ENDDO
	469
[3281]	470	ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') ) THEN
	471	DO i = nxlg, nxrg
	472	DO j = nysg, nyng
	473	chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
	474	ENDDO
[3228]	475	ENDDO
[3281]	476	CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )
	477	ENDIF
	478	ENDDO
	479
	480	ENDIF
[3228]	481
[3281]	482	END SUBROUTINE chem_3d_data_averaging
	483
	484	!
[3228]	485	!------------------------------------------------------------------------------!
	486	!
	487	! Description:
	488	! ------------
[2535]	489	!> Subroutine to initialize and set all boundary conditions for chemical species
[2425]	490	!------------------------------------------------------------------------------!
	491
[3281]	492	SUBROUTINE chem_boundary_conds( mode )
	493
	494	USE control_parameters, &
	495	ONLY: air_chemistry, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
	496	bc_radiation_s
	497	USE indices, &
	498	ONLY: nxl, nxr, nxlg, nxrg, nyng, nysg, nzt
	499
[2425]	500
[3281]	501	USE arrays_3d, &
	502	ONLY: dzu
	503	USE surface_mod, &
	504	ONLY: bc_h
[2425]	505
[3281]	506	CHARACTER (len=*), INTENT(IN) :: mode
	507	INTEGER(iwp) :: i !< grid index x direction.
	508	INTEGER(iwp) :: j !< grid index y direction.
	509	INTEGER(iwp) :: k !< grid index z direction.
	510	INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing.
	511	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls.
	512	INTEGER(iwp) :: m !< running index surface elements.
	513	INTEGER(iwp) :: lsp !< running index for chem spcs.
	514	INTEGER(iwp) :: lph !< running index for photolysis frequencies
[2425]	515
[2615]	516
[3281]	517	SELECT CASE ( TRIM( mode ) )
	518	CASE ( 'init' )
[2626]	519
[3281]	520	IF ( bc_cs_b == 'dirichlet' ) THEN
	521	ibc_cs_b = 0
	522	ELSEIF ( bc_cs_b == 'neumann' ) THEN
	523	ibc_cs_b = 1
	524	ELSE
	525	message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'
	526	CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 ) !<
	527	ENDIF
[2425]	528	!
[3281]	529	!-- Set Integer flags and check for possible erroneous settings for top
	530	!-- boundary condition.
	531	IF ( bc_cs_t == 'dirichlet' ) THEN
	532	ibc_cs_t = 0
	533	ELSEIF ( bc_cs_t == 'neumann' ) THEN
	534	ibc_cs_t = 1
	535	ELSEIF ( bc_cs_t == 'initial_gradient' ) THEN
	536	ibc_cs_t = 2
	537	ELSEIF ( bc_cs_t == 'nested' ) THEN
	538	ibc_cs_t = 3
	539	ELSE
	540	message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'
	541	CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
	542	ENDIF
[2425]	543
[3281]	544
	545	CASE ( 'set_bc_bottomtop' )
	546	!-- Bottom boundary condtions for chemical species
	547	DO lsp = 1, nspec
	548	IF ( ibc_cs_b == 0 ) THEN
	549	DO l = 0, 1
	550	!-- Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
[3282]	551	!-- the chem_species(nspec)%conc_p value at the topography top (k-1)
[3281]	552	!-- is set; for downward-facing surfaces (l=1), kb=1, i.e. the
	553	!-- value at the topography bottom (k+1) is set.
[2425]	554
[3281]	555	kb = MERGE( -1, 1, l == 0 )
	556	!$OMP PARALLEL DO PRIVATE( i, j, k )
	557	DO m = 1, bc_h(l)%ns
	558	i = bc_h(l)%i(m)
	559	j = bc_h(l)%j(m)
	560	k = bc_h(l)%k(m)
	561	chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i)
	562	ENDDO
	563	ENDDO
	564
	565	ELSEIF ( ibc_cs_b == 1 ) THEN
	566	!-- in boundary_conds there is som extra loop over m here for passive tracer
	567	DO l = 0, 1
	568	kb = MERGE( -1, 1, l == 0 )
	569	!$OMP PARALLEL DO PRIVATE( i, j, k )
	570	DO m = 1, bc_h(l)%ns
	571	i = bc_h(l)%i(m)
	572	j = bc_h(l)%j(m)
	573	k = bc_h(l)%k(m)
	574	chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i)
[2535]	575
[3281]	576	ENDDO
[2615]	577	ENDDO
[3281]	578	ENDIF
[3287]	579	ENDDO ! end lsp loop
[3281]	580
[3287]	581	!-- Top boundary conditions for chemical species - Should this not be done for all species?
[3281]	582	IF ( ibc_cs_t == 0 ) THEN
	583	DO lsp = 1, nspec
	584	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)
[2615]	585	ENDDO
[3281]	586	ELSEIF ( ibc_cs_t == 1 ) THEN
	587	DO lsp = 1, nspec
	588	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
	589	ENDDO
	590	ELSEIF ( ibc_cs_t == 2 ) THEN
	591	DO lsp = 1, nspec
	592	chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
	593	ENDDO
[2615]	594	ENDIF
[2425]	595	!
[3281]	596	CASE ( 'set_bc_lateral' )
	597	!-- Lateral boundary conditions for chem species at inflow boundary
	598	!-- are automatically set when chem_species concentration is
	599	!-- initialized. The initially set value at the inflow boundary is not
	600	!-- touched during time integration, hence, this boundary value remains
	601	!-- at a constant value, which is the concentration that flows into the
	602	!-- domain.
	603	!-- Lateral boundary conditions for chem species at outflow boundary
[2535]	604
[3281]	605	IF ( bc_radiation_s ) THEN
	606	DO lsp = 1, nspec
	607	chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
	608	ENDDO
	609	ELSEIF ( bc_radiation_n ) THEN
	610	DO lsp = 1, nspec
	611	chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
	612	ENDDO
	613	ELSEIF ( bc_radiation_l ) THEN
	614	DO lsp = 1, nspec
	615	chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
	616	ENDDO
	617	ELSEIF ( bc_radiation_r ) THEN
	618	DO lsp = 1, nspec
	619	chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
	620	ENDDO
	621	ENDIF
	622
	623	END SELECT
[2535]	624
[3281]	625	END SUBROUTINE chem_boundary_conds
[2831]	626
[2425]	627	!
[2535]	628	!------------------------------------------------------------------------------!
[2831]	629	! Description:
	630	! ------------
	631	!> Boundary conditions for prognostic variables in chemistry: decycling in the
	632	!> x-direction
	633	!-----------------------------------------------------------------------------
[3281]	634	SUBROUTINE chem_boundary_conds_decycle (cs_3d, cs_pr_init )
	635	USE pegrid, &
	636	ONLY: myid
[2831]	637
[3281]	638	IMPLICIT NONE
	639	INTEGER(iwp) :: boundary !<
	640	INTEGER(iwp) :: ee !<
	641	INTEGER(iwp) :: copied !<
	642	INTEGER(iwp) :: i !<
	643	INTEGER(iwp) :: j !<
	644	INTEGER(iwp) :: k !<
	645	INTEGER(iwp) :: ss !<
	646	REAL(wp), DIMENSION(nzb:nzt+1) :: cs_pr_init
	647	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: cs_3d
	648	REAL(wp) :: flag !< flag to mask topography grid points
[2831]	649
[3281]	650	flag = 0.0_wp
[2854]	651
[3281]	652
	653	!-- Left and right boundaries
	654	IF ( decycle_chem_lr .AND. bc_lr_cyc ) THEN
[2854]	655
[3281]	656	DO boundary = 1, 2
[2854]	657
[3281]	658	IF ( decycle_method(boundary) == 'dirichlet' ) THEN
[2854]	659	!
[3281]	660	!-- Initial profile is copied to ghost and first three layers
	661	ss = 1
	662	ee = 0
	663	IF ( boundary == 1 .AND. nxl == 0 ) THEN
	664	ss = nxlg
	665	ee = nxl+2
	666	ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN
	667	ss = nxr-2
	668	ee = nxrg
	669	ENDIF
[2854]	670
[3281]	671	DO i = ss, ee
	672	DO j = nysg, nyng
	673	DO k = nzb+1, nzt
	674	flag = MERGE( 1.0_wp, 0.0_wp, &
	675	BTEST( wall_flags_0(k,j,i), 0 ) )
	676	cs_3d(k,j,i) = cs_pr_init(k) * flag
	677	ENDDO
[2854]	678	ENDDO
[2831]	679	ENDDO
	680
[3281]	681	ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN
[2854]	682	!
[3281]	683	!-- The value at the boundary is copied to the ghost layers to simulate
	684	!-- an outlet with zero gradient
	685	ss = 1
	686	ee = 0
	687	IF ( boundary == 1 .AND. nxl == 0 ) THEN
	688	ss = nxlg
	689	ee = nxl-1
	690	copied = nxl
	691	ELSEIF ( boundary == 2 .AND. nxr == nx ) THEN
	692	ss = nxr+1
	693	ee = nxrg
	694	copied = nxr
	695	ENDIF
[2854]	696
[3281]	697	DO i = ss, ee
	698	DO j = nysg, nyng
	699	DO k = nzb+1, nzt
	700	flag = MERGE( 1.0_wp, 0.0_wp, &
	701	BTEST( wall_flags_0(k,j,i), 0 ) )
	702	cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
	703	ENDDO
[2854]	704	ENDDO
[2831]	705	ENDDO
[2854]	706
[3281]	707	ELSE
	708	WRITE(message_string,*) &
	709	'unknown decycling method: decycle_method (', &
	710	boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
	711	CALL message( 'chem_boundary_conds_decycle', 'CM0431', &
	712	1, 2, 0, 6, 0 )
	713	ENDIF
	714	ENDDO
	715	ENDIF
[2831]	716
[3281]	717
	718	!-- South and north boundaries
	719	IF ( decycle_chem_ns .AND. bc_ns_cyc ) THEN
[2854]	720
[3281]	721	DO boundary = 3, 4
[2854]	722
[3281]	723	IF ( decycle_method(boundary) == 'dirichlet' ) THEN
[2854]	724	!
[3281]	725	!-- Initial profile is copied to ghost and first three layers
	726	ss = 1
	727	ee = 0
	728	IF ( boundary == 3 .AND. nys == 0 ) THEN
	729	ss = nysg
	730	ee = nys+2
	731	ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN
	732	ss = nyn-2
	733	ee = nyng
	734	ENDIF
[2854]	735
[3281]	736	DO i = nxlg, nxrg
	737	DO j = ss, ee
	738	DO k = nzb+1, nzt
	739	flag = MERGE( 1.0_wp, 0.0_wp, &
	740	BTEST( wall_flags_0(k,j,i), 0 ) )
	741	cs_3d(k,j,i) = cs_pr_init(k) * flag
	742	ENDDO
[2854]	743	ENDDO
[2831]	744	ENDDO
[3281]	745
	746
	747	ELSEIF ( decycle_method(boundary) == 'neumann' ) THEN
[2854]	748	!
[3281]	749	!-- The value at the boundary is copied to the ghost layers to simulate
	750	!-- an outlet with zero gradient
	751	ss = 1
	752	ee = 0
	753	IF ( boundary == 3 .AND. nys == 0 ) THEN
	754	ss = nysg
	755	ee = nys-1
	756	copied = nys
	757	ELSEIF ( boundary == 4 .AND. nyn == ny ) THEN
	758	ss = nyn+1
	759	ee = nyng
	760	copied = nyn
	761	ENDIF
[2831]	762
[3281]	763	DO i = nxlg, nxrg
	764	DO j = ss, ee
	765	DO k = nzb+1, nzt
	766	flag = MERGE( 1.0_wp, 0.0_wp, &
	767	BTEST( wall_flags_0(k,j,i), 0 ) )
	768	cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
	769	ENDDO
[2854]	770	ENDDO
[2831]	771	ENDDO
[2854]	772
[3281]	773	ELSE
	774	WRITE(message_string,*) &
	775	'unknown decycling method: decycle_method (', &
	776	boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
	777	CALL message( 'chem_boundary_conds_decycle', 'CM0432', &
	778	1, 2, 0, 6, 0 )
	779	ENDIF
	780	ENDDO
	781	ENDIF
	782	END SUBROUTINE chem_boundary_conds_decycle
	783	!
[3085]	784	!------------------------------------------------------------------------------!
[2831]	785	!
[3085]	786	! Description:
	787	! ------------
	788	!> Subroutine for checking data output for chemical species
[2831]	789	!------------------------------------------------------------------------------!
[3085]	790
[3281]	791	SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
[3085]	792
	793
[3281]	794	USE control_parameters, &
	795	ONLY: data_output, message_string
[3085]	796
[3281]	797	IMPLICIT NONE
[3085]	798
[3281]	799	CHARACTER (LEN=*) :: unit !<
	800	CHARACTER (LEN=*) :: var !<
[3085]	801
[3281]	802	INTEGER(iwp) :: i, lsp
	803	INTEGER(iwp) :: ilen
	804	INTEGER(iwp) :: k
[3085]	805
[3281]	806	CHARACTER(len=16) :: spec_name
[3085]	807
[3281]	808	unit = 'illegal'
[3085]	809
[3281]	810	spec_name = TRIM(var(4:)) !< var 1:3 is 'kc_'.
[3085]	811
[3281]	812	DO lsp=1,nspec
	813	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN
	814	unit = 'ppm'
	815	ENDIF
	816	! It is possible to plant PM10 and PM25 into the gasphase chemistry code
	817	! as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
	818	! set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
	819	IF (spec_name(1:2) == 'PM') THEN
[3276]	820	unit = 'kg m-3'
[3281]	821	ENDIF
	822	ENDDO
[3085]	823
[3281]	824	DO lsp=1,nphot
	825	IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name)) THEN
	826	unit = 'sec-1'
	827	ENDIF
	828	ENDDO
[3085]	829
	830
	831	RETURN
[3281]	832	END SUBROUTINE chem_check_data_output
[2831]	833	!
[3085]	834	!------------------------------------------------------------------------------!
	835	!
[2535]	836	! Description:
	837	! ------------
[3085]	838	!> Subroutine for checking data output of profiles for chemistry model
[2535]	839	!------------------------------------------------------------------------------!
[2615]	840
[3281]	841	SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
[3085]	842
[3281]	843	USE arrays_3d
	844	USE control_parameters, &
	845	ONLY: data_output_pr, message_string, air_chemistry
	846	USE indices
	847	USE profil_parameter
	848	USE statistics
[3085]	849
	850
[3281]	851	IMPLICIT NONE
[3085]	852
[3281]	853	CHARACTER (LEN=*) :: unit !<
	854	CHARACTER (LEN=*) :: variable !<
	855	CHARACTER (LEN=*) :: dopr_unit
	856	CHARACTER(len=16) :: spec_name
[3228]	857
[3281]	858	INTEGER(iwp) :: var_count, lsp !<
	859
[3085]	860
[3281]	861	spec_name = TRIM(variable(4:))
[3085]	862
	863	IF ( .NOT. air_chemistry ) THEN
[3281]	864	message_string = 'data_output_pr = ' // &
	865	TRIM( data_output_pr(var_count) ) // ' is not imp' // &
	866	'lemented for air_chemistry = .FALSE.'
	867	CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )
[3085]	868
	869	ELSE
	870	DO lsp = 1, nspec
	871	IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
[3281]	872	cs_pr_count = cs_pr_count+1
	873	cs_pr_index(cs_pr_count) = lsp
	874	dopr_index(var_count) = pr_palm+cs_pr_count
	875	dopr_unit = 'ppm'
	876	IF (spec_name(1:2) == 'PM') THEN
[3276]	877	dopr_unit = 'kg m-3'
[3281]	878	ENDIF
[3085]	879	hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
	880	unit = dopr_unit
	881	ENDIF
	882	ENDDO
	883	ENDIF
	884
[3281]	885	END SUBROUTINE chem_check_data_output_pr
[3085]	886
[3281]	887	!
	888	!------------------------------------------------------------------------------!
	889	! Description:
	890	! ------------
	891	!> Check parameters routine for chemistry_model_mod
	892	!------------------------------------------------------------------------------!
	893	SUBROUTINE chem_check_parameters
[3085]	894
[3281]	895	IMPLICIT NONE
	896
	897	LOGICAL :: found
	898	INTEGER (iwp) :: lsp_usr !< running index for user defined chem spcs
	899	INTEGER (iwp) :: lsp !< running index for chem spcs.
	900
	901
	902	!!-- check for chemical reactions status
	903	IF ( chem_gasphase_on ) THEN
	904	message_string = 'Chemical reactions: ON'
	905	CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
	906	ELSEIF ( .not. (chem_gasphase_on) ) THEN
	907	message_string = 'Chemical reactions: OFF'
	908	CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
	909	ENDIF
	910
	911	!-- check for chemistry time-step
	912	IF ( call_chem_at_all_substeps ) THEN
	913	message_string = 'Chemistry is calculated at all meteorology time-step'
	914	CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
	915	ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN
	916	message_string = 'Sub-time-steps are skipped for chemistry time-steps'
	917	CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
	918	ENDIF
	919
	920	!-- check for photolysis scheme
	921	IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant') ) THEN
	922	message_string = 'Incorrect photolysis scheme selected, please check spelling'
	923	CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
	924	ENDIF
	925
	926	!-- check for decycling of chem species
	927	IF ( (.not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) ) THEN
[3318]	928	message_string = 'Incorrect boundary conditions. Only neumann or ' &
	929	// 'dirichlet &available for decycling chemical species '
[3281]	930	CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
	931	ENDIF
	932
	933	!-- check for initial chem species input
	934	lsp_usr = 1
	935	lsp = 1
	936	DO WHILE ( cs_name (lsp_usr) /= 'novalue')
	937	found = .FALSE.
	938	DO lsp = 1, nvar
	939	IF ( TRIM(cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
	940	found = .TRUE.
	941	EXIT
	942	ENDIF
	943	ENDDO
	944	IF ( .not. found ) THEN
	945	message_string = 'Incorrect input for initial surface vaue: ' // TRIM(cs_name(lsp_usr))
	946	CALL message( 'chem_check_parameters', 'CM0427', 0, 1, 0, 6, 0 )
	947	ENDIF
	948	lsp_usr = lsp_usr + 1
	949	ENDDO
	950
	951	!-- check for surface emission flux chem species
	952
	953	lsp_usr = 1
	954	lsp = 1
	955	DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
	956	found = .FALSE.
	957	DO lsp = 1, nvar
	958	IF ( TRIM(surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
	959	found = .TRUE.
	960	EXIT
	961	ENDIF
	962	ENDDO
	963	IF ( .not. found ) THEN
	964	message_string = 'Incorrect input of chemical species for surface emission fluxes: ' &
	965	// TRIM(surface_csflux_name(lsp_usr))
	966	CALL message( 'chem_check_parameters', 'CM0428', 0, 1, 0, 6, 0 )
	967	ENDIF
	968	lsp_usr = lsp_usr + 1
	969	ENDDO
	970
	971	END SUBROUTINE chem_check_parameters
	972
[3085]	973	!
	974	!------------------------------------------------------------------------------!
	975	!
	976	! Description:
	977	! ------------
	978	!> Subroutine defining 2D output variables for chemical species
[3298]	979	!> @todo: Remove "mode" from argument list, not used.
[3085]	980	!------------------------------------------------------------------------------!
	981
[3281]	982	SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf, &
	983	two_d, nzb_do, nzt_do, fill_value )
	984
	985	USE indices
[3085]	986
[3281]	987	USE kinds
[3085]	988
[3281]	989	USE pegrid, ONLY: myid, threads_per_task
[3085]	990
[3281]	991	IMPLICIT NONE
[3085]	992
[3281]	993	CHARACTER (LEN=*) :: grid !<
	994	CHARACTER (LEN=*) :: mode !<
	995	CHARACTER (LEN=*) :: variable !<
[3282]	996	INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data
[3281]	997	INTEGER(iwp) :: nzb_do !< lower limit of the domain (usually nzb)
	998	INTEGER(iwp) :: nzt_do !< upper limit of the domain (usually nzt+1)
	999	LOGICAL :: found !<
	1000	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
	1001	REAL(wp) :: fill_value
	1002	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) :: local_pf !<
[3085]	1003
[3281]	1004	!-- local variables.
	1005	CHARACTER(len=16) :: spec_name
	1006	INTEGER(iwp) :: lsp
[3282]	1007	INTEGER(iwp) :: i !< grid index along x-direction
	1008	INTEGER(iwp) :: j !< grid index along y-direction
	1009	INTEGER(iwp) :: k !< grid index along z-direction
	1010	INTEGER(iwp) :: m !< running index surface elements
	1011	INTEGER(iwp) :: char_len !< length of a character string
[3281]	1012	found = .TRUE.
	1013	char_len = LEN_TRIM(variable)
[3085]	1014
[3281]	1015	spec_name = TRIM( variable(4:char_len-3) )
[3085]	1016
[3281]	1017	DO lsp=1,nspec
	1018	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) .AND. &
	1019	( (variable(char_len-2:) == '_xy') .OR. &
	1020	(variable(char_len-2:) == '_xz') .OR. &
	1021	(variable(char_len-2:) == '_yz') ) ) THEN
[3373]	1022	!
	1023	!-- todo: remove or replace by "CALL message" mechanism (kanani)
	1024	! IF(myid == 0) WRITE(6,*) 'Output of species ' // TRIM(variable) // &
	1025	! TRIM(chem_species(lsp)%name)
[3281]	1026	IF (av == 0) THEN
	1027	DO i = nxl, nxr
	1028	DO j = nys, nyn
	1029	DO k = nzb_do, nzt_do
	1030	local_pf(i,j,k) = MERGE( &
	1031	chem_species(lsp)%conc(k,j,i), &
	1032	REAL( fill_value, KIND = wp ), &
	1033	BTEST( wall_flags_0(k,j,i), 0 ) )
[3085]	1034
	1035
[3281]	1036	ENDDO
[3085]	1037	ENDDO
	1038	ENDDO
[3281]	1039
	1040	ELSE
	1041	DO i = nxl, nxr
	1042	DO j = nys, nyn
	1043	DO k = nzb_do, nzt_do
	1044	local_pf(i,j,k) = MERGE( &
	1045	chem_species(lsp)%conc(k,j,i), &
	1046	REAL( fill_value, KIND = wp ), &
	1047	BTEST( wall_flags_0(k,j,i), 0 ) )
	1048	ENDDO
[3085]	1049	ENDDO
	1050	ENDDO
[3281]	1051	ENDIF
	1052	grid = 'zu'
[2535]	1053	ENDIF
[3281]	1054	ENDDO
[3085]	1055
[3281]	1056	RETURN
	1057
	1058	END SUBROUTINE chem_data_output_2d
[3085]	1059
	1060	!
	1061	!------------------------------------------------------------------------------!
	1062	!
	1063	! Description:
	1064	! ------------
	1065	!> Subroutine defining 3D output variables for chemical species
	1066	!------------------------------------------------------------------------------!
	1067
[3281]	1068	SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
[3085]	1069
	1070
[3281]	1071	USE indices
[3228]	1072
[3281]	1073	USE kinds
[3085]	1074
	1075
[3281]	1076	IMPLICIT NONE
[3085]	1077
[3282]	1078	CHARACTER (LEN=*) :: variable !<
	1079	INTEGER(iwp) :: av !<
	1080	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
	1081	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
[3085]	1082
[3282]	1083	LOGICAL :: found !<
[3085]	1084
[3281]	1085	REAL(wp) :: fill_value !<
	1086	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf
[3228]	1087
[3085]	1088
[3281]	1089	!-- local variables
[3085]	1090
[3281]	1091	INTEGER :: i, j, k, lsp
	1092	CHARACTER(len=16) :: spec_name
[3085]	1093
	1094
[3281]	1095	found = .FALSE.
[3085]	1096
[3281]	1097	spec_name = TRIM(variable(4:))
[3085]	1098
	1099	DO lsp=1,nspec
[3114]	1100	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)) THEN
[3373]	1101	!
	1102	!-- todo: remove or replace by "CALL message" mechanism (kanani)
	1103	! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM(variable) // &
	1104	! TRIM(chem_species(lsp)%name)
[3228]	1105
[3085]	1106	IF (av == 0) THEN
	1107	DO i = nxl, nxr
	1108	DO j = nys, nyn
[3114]	1109	DO k = nzb_do, nzt_do
[3085]	1110	local_pf(i,j,k) = MERGE( &
	1111	chem_species(lsp)%conc(k,j,i), &
	1112	REAL( fill_value, KIND = wp ), &
	1113	BTEST( wall_flags_0(k,j,i), 0 ) )
	1114	ENDDO
	1115	ENDDO
	1116	ENDDO
	1117
	1118	ELSE
	1119	DO i = nxl, nxr
	1120	DO j = nys, nyn
[3114]	1121	DO k = nzb_do, nzt_do
[3085]	1122	local_pf(i,j,k) = MERGE( &
	1123	chem_species(lsp)%conc_av(k,j,i),&
	1124	REAL( fill_value, KIND = wp ), &
	1125	BTEST( wall_flags_0(k,j,i), 0 ) )
	1126	ENDDO
	1127	ENDDO
	1128	ENDDO
	1129	ENDIF
	1130
[3281]	1131	found = .TRUE.
[3085]	1132	ENDIF
	1133	ENDDO
	1134
	1135	RETURN
[3281]	1136	END SUBROUTINE chem_data_output_3d
[3085]	1137	!
	1138	!------------------------------------------------------------------------------!
	1139	!
	1140	! Description:
	1141	! ------------
	1142	!> Subroutine defining mask output variables for chemical species
	1143	!------------------------------------------------------------------------------!
	1144
[3281]	1145	SUBROUTINE chem_data_output_mask( av, variable, found, local_pf )
	1146
	1147	USE control_parameters
	1148	USE indices
	1149	USE kinds
	1150	USE pegrid, ONLY: myid, threads_per_task
[3085]	1151
	1152
[3281]	1153	IMPLICIT NONE
[3085]	1154
[3282]	1155	CHARACTER (LEN=*):: variable !<
	1156	INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data
	1157	LOGICAL :: found !<
[3281]	1158	REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
	1159	local_pf !<
[3085]	1160
	1161
[3281]	1162	!-- local variables.
	1163	CHARACTER(len=16) :: spec_name
	1164	INTEGER(iwp) :: lsp
[3282]	1165	INTEGER(iwp) :: i !< grid index along x-direction
	1166	INTEGER(iwp) :: j !< grid index along y-direction
	1167	INTEGER(iwp) :: k !< grid index along z-direction
[3281]	1168	found = .TRUE.
[3085]	1169
[3281]	1170	spec_name = TRIM( variable(4:) )
	1171	!av == 0
[3085]	1172
	1173	DO lsp=1,nspec
	1174	IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) ) THEN
[3373]	1175	!
	1176	!-- todo: remove or replace by "CALL message" mechanism (kanani)
	1177	! IF(myid == 0 .AND. chem_debug0 ) WRITE(6,*) 'Output of species ' // TRIM(variable) // &
	1178	! TRIM(chem_species(lsp)%name)
[3085]	1179	IF (av == 0) THEN
	1180	DO i = 1, mask_size_l(mid,1)
	1181	DO j = 1, mask_size_l(mid,2)
	1182	DO k = 1, mask_size(mid,3)
	1183	local_pf(i,j,k) = chem_species(lsp)%conc(mask_k(mid,k), &
	1184	mask_j(mid,j), mask_i(mid,i))
	1185	ENDDO
	1186	ENDDO
	1187	ENDDO
	1188
	1189	ELSE
	1190	DO i = 1, mask_size_l(mid,1)
	1191	DO j = 1, mask_size_l(mid,2)
	1192	DO k = 1, mask_size_l(mid,3)
	1193	local_pf(i,j,k) = chem_species(lsp)%conc_av(mask_k(mid,k), &
	1194	mask_j(mid,j), mask_i(mid,i))
	1195	ENDDO
	1196	ENDDO
	1197	ENDDO
	1198	ENDIF
	1199	found = .FALSE.
	1200	ENDIF
	1201	ENDDO
	1202
	1203	RETURN
[3281]	1204
	1205	END SUBROUTINE chem_data_output_mask
[3085]	1206
	1207	!
	1208	!------------------------------------------------------------------------------!
	1209	!
	1210	! Description:
	1211	! ------------
	1212	!> Subroutine defining appropriate grid for netcdf variables.
	1213	!> It is called out from subroutine netcdf.
	1214	!------------------------------------------------------------------------------!
[3281]	1215	SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
[3085]	1216
[3281]	1217	IMPLICIT NONE
[3085]	1218
[3282]	1219	CHARACTER (LEN=*), INTENT(IN) :: var !<
	1220	LOGICAL, INTENT(OUT) :: found !<
	1221	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
	1222	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
	1223	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
[3085]	1224
[3281]	1225	found = .TRUE.
[3085]	1226
[3281]	1227	IF (var(1:3) == 'kc_') THEN !< always the same grid for chemistry variables
	1228	grid_x = 'x'
	1229	grid_y = 'y'
	1230	grid_z = 'zu'
	1231	ELSE
	1232	found = .FALSE.
	1233	grid_x = 'none'
	1234	grid_y = 'none'
	1235	grid_z = 'none'
	1236	ENDIF
[3085]	1237
	1238
[3281]	1239	END SUBROUTINE chem_define_netcdf_grid
[3085]	1240	!
	1241	!------------------------------------------------------------------------------!
	1242	!
	1243	! Description:
	1244	! ------------
	1245	!> Subroutine defining header output for chemistry model
	1246	!------------------------------------------------------------------------------!
[3281]	1247	SUBROUTINE chem_header ( io )
[3085]	1248
[3281]	1249	IMPLICIT NONE
[3085]	1250
	1251	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
[3281]	1252	INTEGER(iwp) :: lsp !< running index for chem spcs
[3282]	1253	INTEGER(iwp) :: cs_fixed
	1254	CHARACTER (LEN=80) :: docsflux_chr
	1255	CHARACTER (LEN=80) :: docsinit_chr
	1256
[3281]	1257	!
	1258	!-- Write chemistry model header
	1259	WRITE( io, 1 )
[3085]	1260
[3281]	1261	!-- Gasphase reaction status
	1262	IF ( chem_gasphase_on ) THEN
	1263	WRITE( io, 2 )
	1264	ELSE
	1265	WRITE( io, 3 )
	1266	ENDIF
	1267
	1268	! Chemistry time-step
	1269	WRITE ( io, 4 ) cs_time_step
	1270
	1271	!-- Emission mode info
	1272	IF ( mode_emis == "DEFAULT" ) THEN
	1273	WRITE( io, 5 )
	1274	ELSEIF ( mode_emis == "PARAMETERIZED" ) THEN
	1275	WRITE( io, 6 )
	1276	ELSEIF ( mode_emis == "PRE-PROCESSED" ) THEN
	1277	WRITE( io, 7 )
	1278	ENDIF
	1279
	1280	!-- Photolysis scheme info
	1281	IF ( photolysis_scheme == "simple" ) THEN
	1282	WRITE( io, 8 )
	1283	ELSEIF (photolysis_scheme == "conastant" ) THEN
	1284	WRITE( io, 9 )
	1285	ENDIF
	1286
	1287	!-- Emission flux info
	1288	lsp = 1
	1289	docsflux_chr ='Chemical species for surface emission flux: '
	1290	DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
	1291	docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ','
	1292	IF ( LEN_TRIM( docsflux_chr ) >= 75 ) THEN
	1293	WRITE ( io, 10 ) docsflux_chr
	1294	docsflux_chr = ' '
	1295	ENDIF
	1296	lsp = lsp + 1
	1297	ENDDO
	1298
	1299	IF ( docsflux_chr /= '' ) THEN
	1300	WRITE ( io, 10 ) docsflux_chr
	1301	ENDIF
	1302
	1303
	1304	!-- initializatoin of Surface and profile chemical species
	1305
	1306	lsp = 1
	1307	docsinit_chr ='Chemical species for initial surface and profile emissions: '
	1308	DO WHILE ( cs_name(lsp) /= 'novalue' )
	1309	docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ','
	1310	IF ( LEN_TRIM( docsinit_chr ) >= 75 ) THEN
	1311	WRITE ( io, 11 ) docsinit_chr
	1312	docsinit_chr = ' '
	1313	ENDIF
	1314	lsp = lsp + 1
	1315	ENDDO
	1316
	1317	IF ( docsinit_chr /= '' ) THEN
	1318	WRITE ( io, 11 ) docsinit_chr
	1319	ENDIF
	1320
	1321	!-- number of variable and fix chemical species and number of reactions
[3282]	1322	cs_fixed = nspec - nvar
	1323	WRITE ( io, * ) ' --> Chemical species, variable: ', nvar
	1324	WRITE ( io, * ) ' --> Chemical species, fixed : ', cs_fixed
	1325	WRITE ( io, * ) ' --> Total number of reactions : ', nreact
[3281]	1326
	1327
	1328	1 FORMAT (//' Chemistry model information:'/ &
	1329	' ----------------------------'/)
	1330	2 FORMAT (' --> Chemical reactions are turned on')
	1331	3 FORMAT (' --> Chemical reactions are turned off')
	1332	4 FORMAT (' --> Time-step for chemical species: ',F6.2, ' s')
	1333	5 FORMAT (' --> Emission mode = DEFAULT ')
	1334	6 FORMAT (' --> Emission mode = PARAMETERIZED ')
	1335	7 FORMAT (' --> Emission mode = PRE-PROCESSED ')
	1336	8 FORMAT (' --> Photolysis scheme used = simple ')
	1337	9 FORMAT (' --> Photolysis scheme used = constant ')
[3282]	1338	10 FORMAT (/' ',A)
	1339	11 FORMAT (/' ',A)
[3085]	1340	!
	1341	!
[3281]	1342	END SUBROUTINE chem_header
[3085]	1343
	1344	!
	1345	!------------------------------------------------------------------------------!
	1346	!
	1347	! Description:
	1348	! ------------
[2535]	1349	!> Subroutine initializating chemistry_model_mod
[2425]	1350	!------------------------------------------------------------------------------!
[3281]	1351	SUBROUTINE chem_init
[2535]	1352
[3281]	1353	USE control_parameters, &
	1354	ONLY: message_string, io_blocks, io_group, turbulent_inflow
	1355	USE arrays_3d, &
	1356	ONLY: mean_inflow_profiles
	1357	USE pegrid
[3085]	1358
[3281]	1359	IMPLICIT none
	1360	!-- local variables
	1361	INTEGER :: i,j !< running index for for horiz numerical grid points
	1362	INTEGER :: lsp !< running index for chem spcs
	1363	INTEGER :: lpr_lev !< running index for chem spcs profile level
[2535]	1364	!
	1365	!-- NOPOINTER version not implemented yet
[2828]	1366	! #if defined( __nopointer )
	1367	! message_string = 'The chemistry module only runs with POINTER version'
	1368	! CALL message( 'chemistry_model_mod', 'CM0001', 1, 2, 0, 6, 0 )
	1369	! #endif
	1370	!
	1371	!-- Allocate memory for chemical species
[3281]	1372	ALLOCATE( chem_species(nspec) )
	1373	ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
	1374	ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
	1375	ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
	1376	ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
	1377	ALLOCATE( phot_frequen(nphot) )
	1378	ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
	1379	ALLOCATE( bc_cs_t_val(nspec) )
[2828]	1380	!
	1381	!-- Initialize arrays
[3281]	1382	spec_conc_1 (:,:,:,:) = 0.0_wp
	1383	spec_conc_2 (:,:,:,:) = 0.0_wp
	1384	spec_conc_3 (:,:,:,:) = 0.0_wp
	1385	spec_conc_av(:,:,:,:) = 0.0_wp
[2535]	1386
[2828]	1387
[3281]	1388	DO lsp = 1, nspec
	1389	chem_species(lsp)%name = spc_names(lsp)
[2535]	1390
[3281]	1391	chem_species(lsp)%conc (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1 (:,:,:,lsp)
	1392	chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2 (:,:,:,lsp)
	1393	chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_3 (:,:,:,lsp)
	1394	chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_av(:,:,:,lsp)
[2425]	1395
[3281]	1396	ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))
	1397	chem_species(lsp)%cssws_av = 0.0_wp
[2828]	1398	!
[3282]	1399	!-- The following block can be useful when emission module is not applied. &
	1400	!-- if emission module is applied the following block will be overwritten.
	1401	ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))
	1402	ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))
	1403	ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
	1404	ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
	1405	chem_species(lsp)%flux_s_cs = 0.0_wp
	1406	chem_species(lsp)%flux_l_cs = 0.0_wp
	1407	chem_species(lsp)%diss_s_cs = 0.0_wp
	1408	chem_species(lsp)%diss_l_cs = 0.0_wp
[2828]	1409	!
	1410	!-- Allocate memory for initial concentration profiles
	1411	!-- (concentration values come from namelist)
[3282]	1412	!-- (@todo (FK): Because of this, chem_init is called in palm before
[2828]	1413	!-- check_parameters, since conc_pr_init is used there.
	1414	!-- We have to find another solution since chem_init should
	1415	!-- eventually be called from init_3d_model!!)
[3281]	1416	ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
	1417	chem_species(lsp)%conc_pr_init(:) = 0.0_wp
[2425]	1418
[3281]	1419	ENDDO
[2425]	1420
[2615]	1421	!
[3281]	1422	!-- Initial concentration of profiles is prescribed by parameters cs_profile
	1423	!-- and cs_heights in the namelist &chemistry_parameters
	1424	CALL chem_init_profiles
	1425
	1426
[2615]	1427	!
	1428	!-- Initialize model variables
[2425]	1429
[3114]	1430
[3281]	1431	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
	1432	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
[2425]	1433
[3228]	1434
[2615]	1435	!-- First model run of a possible job queue.
[3282]	1436	!-- Initial profiles of the variables must be computed.
[3281]	1437	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
[3282]	1438	CALL location_message( 'initializing with 1D chemistry model profiles', .FALSE. )
[2615]	1439	!
	1440	!-- Transfer initial profiles to the arrays of the 3D model
[3281]	1441	DO lsp = 1, nspec
	1442	DO i = nxlg, nxrg
	1443	DO j = nysg, nyng
	1444	DO lpr_lev = 1, nz + 1
	1445	chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
	1446	ENDDO
[2615]	1447	ENDDO
[3281]	1448	ENDDO
	1449	ENDDO
	1450
	1451	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
	1452	THEN
[3282]	1453	CALL location_message( 'initializing with constant chemistry profiles', .FALSE. )
[2425]	1454
[3281]	1455	DO lsp = 1, nspec
	1456	DO i = nxlg, nxrg
	1457	DO j = nysg, nyng
	1458	chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init
	1459	ENDDO
[2615]	1460	ENDDO
	1461	ENDDO
	1462
[3281]	1463	ENDIF
[2615]	1464
[2535]	1465	!
[3281]	1466	!-- If required, change the surface chem spcs at the start of the 3D run
	1467	IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN
	1468	DO lsp = 1, nspec
	1469	chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) + &
	1470	cs_surface_initial_change(lsp)
	1471	ENDDO
	1472	ENDIF
	1473	!
	1474	!-- Initiale old and new time levels.
	1475	DO lsp = 1, nvar
	1476	chem_species(lsp)%tconc_m = 0.0_wp
	1477	chem_species(lsp)%conc_p = chem_species(lsp)%conc
[2615]	1478	ENDDO
[2592]	1479
[3281]	1480	ENDIF
[2535]	1481
[3114]	1482
	1483
[3281]	1484	!--- new code add above this line
	1485	DO lsp = 1, nphot
	1486	phot_frequen(lsp)%name = phot_names(lsp)
[3373]	1487	!
	1488	!-- todo: remove or replace by "CALL message" mechanism (kanani)
[3281]	1489	! IF( myid == 0 ) THEN
	1490	! WRITE(6,'(a,i4,3x,a)') 'Photolysis: ',lsp,trim(phot_names(lsp))
	1491	! ENDIF
	1492	phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => freq_1(:,:,:,lsp)
	1493	ENDDO
[2425]	1494
[3281]	1495	RETURN
[2425]	1496
[3085]	1497	END SUBROUTINE chem_init
[2425]	1498
[2535]	1499	!
	1500	!------------------------------------------------------------------------------!
	1501	!
	1502	! Description:
	1503	! ------------
[3085]	1504	!> Subroutine defining initial vertical profiles of chemical species (given by
	1505	!> namelist parameters chem_profiles and chem_heights) --> which should work
	1506	!> analogue to parameters u_profile, v_profile and uv_heights)
[2535]	1507	!------------------------------------------------------------------------------!
[3281]	1508	SUBROUTINE chem_init_profiles !< SUBROUTINE is called from chem_init in case of
	1509	!< TRIM( initializing_actions ) /= 'read_restart_data'
	1510	!< We still need to see what has to be done in case of restart run
	1511	USE chem_modules
[2615]	1512
[3281]	1513	IMPLICIT NONE
[2425]	1514
[3281]	1515	!-- Local variables
	1516	INTEGER :: lsp !< running index for number of species in derived data type species_def
	1517	INTEGER :: lsp_usr !< running index for number of species (user defined) in cs_names, cs_profiles etc
	1518	INTEGER :: lpr_lev !< running index for profile level for each chem spcs.
	1519	INTEGER :: npr_lev !< the next available profile lev
[2425]	1520
[3085]	1521	!-----------------
	1522	!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
	1523	!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
	1524	!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
	1525	!-- "cs_heights" are prescribed, their values will!override the constant profile given by
	1526	!-- "cs_surface".
[3114]	1527	!
[3281]	1528	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1529	lsp_usr = 1
	1530	DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' ) !'novalue' is the default
	1531	DO lsp = 1, nspec !
[3282]	1532	!-- create initial profile (conc_pr_init) for each chemical species
[3281]	1533	IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) ) THEN !
	1534	IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN
[3282]	1535	!-- set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
	1536	DO lpr_lev = 0, nzt+1
	1537	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
	1538	ENDDO
[3281]	1539	ELSE
	1540	IF ( cs_heights(1,1) /= 0.0_wp ) THEN
	1541	message_string = 'The surface value of cs_heights must be 0.0'
	1542	CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
	1543	ENDIF
	1544
	1545	use_prescribed_profile_data = .TRUE.
	1546
	1547	npr_lev = 1
	1548	! chem_species(lsp)%conc_pr_init(0) = 0.0_wp
	1549	DO lpr_lev = 1, nz+1
	1550	IF ( npr_lev < 100 ) THEN
	1551	DO WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
	1552	npr_lev = npr_lev + 1
	1553	IF ( npr_lev == 100 ) THEN
	1554	message_string = 'number of chem spcs exceeding the limit'
	1555	CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )
	1556	EXIT
	1557	ENDIF
	1558	ENDDO
	1559	ENDIF
	1560	IF ( npr_lev < 100 .AND. cs_heights(lsp_usr, npr_lev + 1) /= 9999999.9_wp ) THEN
	1561	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) + &
	1562	( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) / &
	1563	( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) * &
	1564	( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
	1565	ELSE
	1566	chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
	1567	ENDIF
	1568	ENDDO
[3085]	1569	ENDIF
[3282]	1570	!-- If a profile is prescribed explicity using cs_profiles and cs_heights, then
	1571	!-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
	1572	!-- on the cs_profiles(lsp_usr,:) and cs_heights(lsp_usr,:).
[3281]	1573	ENDIF
	1574	ENDDO
	1575	lsp_usr = lsp_usr + 1
[3085]	1576	ENDDO
[3281]	1577	ENDIF
[2425]	1578
[3281]	1579	END SUBROUTINE chem_init_profiles
[2828]	1580	!
[2535]	1581	!------------------------------------------------------------------------------!
	1582	!
	1583	! Description:
	1584	! ------------
	1585	!> Subroutine to integrate chemical species in the given chemical mechanism
	1586	!------------------------------------------------------------------------------!
	1587
[3281]	1588	SUBROUTINE chem_integrate_ij (i, j)
[2425]	1589
[3281]	1590	USE cpulog, &
	1591	ONLY: cpu_log, log_point
[3287]	1592	USE statistics, &
[3281]	1593	ONLY: weight_pres
[3287]	1594	USE control_parameters, &
[3281]	1595	ONLY: dt_3d, intermediate_timestep_count,simulated_time
[2425]	1596
[3281]	1597	IMPLICIT none
	1598	INTEGER,INTENT(IN) :: i,j
[2425]	1599
[3281]	1600	!-- local variables
[3282]	1601	INTEGER(iwp) :: lsp !< running index for chem spcs.
	1602	INTEGER(iwp) :: lph !< running index for photolysis frequencies
[3281]	1603	INTEGER :: k,m,istatf
	1604	INTEGER,DIMENSION(20) :: istatus
	1605	REAL(kind=wp),DIMENSION(nzb+1:nzt,nspec) :: tmp_conc
	1606	REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_temp
	1607	REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_qvap
	1608	REAL(kind=wp),DIMENSION(nzb+1:nzt,nphot) :: tmp_phot
	1609	REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_fact
	1610	REAL(kind=wp),DIMENSION(nzb+1:nzt) :: tmp_fact_i !< conversion factor between
	1611	!< molecules cm^{-3} and ppm
[3185]	1612
[3281]	1613	INTEGER,DIMENSION(nzb+1:nzt) :: nacc !< Number of accepted steps
	1614	INTEGER,DIMENSION(nzb+1:nzt) :: nrej !< Number of rejected steps
[3185]	1615
[3281]	1616	REAL(wp) :: conv !< conversion factor
	1617	REAL(wp), PARAMETER :: ppm2fr = 1.0e-6_wp !< Conversion factor ppm to fraction
	1618	REAL(wp), PARAMETER :: xm_air = 28.96_wp !< Mole mass of dry air
	1619	REAL(wp), PARAMETER :: xm_h2o = 18.01528_wp !< Mole mass of water vapor
	1620	REAL(wp), PARAMETER :: pref_i = 1._wp / 100000.0_wp !< inverse reference pressure (1/Pa)
	1621	REAL(wp), PARAMETER :: t_std = 273.15_wp !< standard pressure (Pa)
	1622	REAL(wp), PARAMETER :: p_std = 101325.0_wp !< standard pressure (Pa)
	1623	REAL(wp), PARAMETER :: vmolcm = 22.414e3_wp !< Mole volume (22.414 l) in cm^{-3}
	1624	REAL(wp), PARAMETER :: xna = 6.022e23_wp !< Avogadro number (molecules/mol)
[2425]	1625
[3281]	1626	REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
[2615]	1627
[3185]	1628
[3281]	1629	REAL(kind=wp) :: dt_chem
[2425]	1630
[2535]	1631	CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' )
[3282]	1632	!< set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
[3281]	1633	IF (chem_gasphase_on) THEN
	1634	nacc = 0
	1635	nrej = 0
[2635]	1636
[3287]	1637	tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
[3282]	1638	! ppm to molecules/cm**3
[3287]	1639	! tmp_fact = 1.e-6_wp6.022e23_wp/(22.414_wp1000._wp) * 273.15_wp * hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp )
[3281]	1640	conv = ppm2fr * xna / vmolcm
	1641	tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
	1642	tmp_fact_i = 1.0_wp/tmp_fact
[2615]	1643
[3281]	1644	IF ( humidity ) THEN
[3292]	1645	IF ( bulk_cloud_model ) THEN
[3281]	1646	tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) * xm_air/xm_h2o * tmp_fact(:)
	1647	ELSE
	1648	tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * tmp_fact(:)
	1649	ENDIF
[3114]	1650	ELSE
[3281]	1651	tmp_qvap(:) = 0.01 * tmp_fact(:) !< Constant value for q if water vapor is not computed
[3114]	1652	ENDIF
	1653
[3281]	1654	DO lsp = 1,nspec
	1655	tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
	1656	ENDDO
[2425]	1657
[3281]	1658	DO lph = 1,nphot
	1659	tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)
	1660	ENDDO
[3373]	1661	!
	1662	!-- todo: remove (kanani)
	1663	! IF(myid == 0 .AND. chem_debug0 ) THEN
	1664	! IF (i == 10 .and. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
	1665	! ENDIF
[2425]	1666
[3281]	1667	!-- Compute length of time step
	1668	IF ( call_chem_at_all_substeps ) THEN
	1669	dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
	1670	ELSE
	1671	dt_chem = dt_3d
	1672	ENDIF
[2425]	1673
[3281]	1674	cs_time_step = dt_chem
[2425]	1675
[3281]	1676	CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' )
[2425]	1677
[3281]	1678	IF(maxval(rcntrl) > 0.0) THEN ! Only if rcntrl is set
	1679	IF( simulated_time <= 2*dt_3d) THEN
	1680	rcntrl_local = 0
[3373]	1681	!
	1682	!-- todo: remove (kanani)
	1683	! WRITE(9,'(a,2f10.3)') 'USE Default rcntrl in the first steps ',simulated_time,dt_3d
[3281]	1684	ELSE
	1685	rcntrl_local = rcntrl
	1686	ENDIF
	1687	ELSE
	1688	rcntrl_local = 0
	1689	END IF
[2425]	1690
[3281]	1691	CALL chem_gasphase_integrate (dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
	1692	icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus)
[2425]	1693
[3281]	1694	CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' )
[2425]	1695
[3281]	1696	DO lsp = 1,nspec
	1697	chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
	1698	ENDDO
[3185]	1699
[2535]	1700
[3281]	1701	ENDIF
	1702	CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' )
	1703
[2615]	1704	RETURN
[3281]	1705	END SUBROUTINE chem_integrate_ij
[2535]	1706	!
	1707	!------------------------------------------------------------------------------!
	1708	!
	1709	! Description:
	1710	! ------------
[3085]	1711	!> Subroutine defining parin for &chemistry_parameters for chemistry model
[2535]	1712	!------------------------------------------------------------------------------!
[3281]	1713	SUBROUTINE chem_parin
[3085]	1714
[3281]	1715	USE chem_modules
	1716	USE control_parameters
[3287]	1717
[3281]	1718	USE kinds
	1719	USE pegrid
	1720	USE statistics
	1721
	1722	IMPLICIT NONE
[2425]	1723
[3282]	1724	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
[3281]	1725	CHARACTER (LEN=3) :: cs_prefix
[2425]	1726
[3282]	1727	REAL(wp), DIMENSION(nmaxfixsteps) :: my_steps !< List of fixed timesteps my_step(1) = 0.0 automatic stepping
	1728	INTEGER(iwp) :: i !<
	1729	INTEGER(iwp) :: j !<
	1730	INTEGER(iwp) :: max_pr_cs_tmp !<
[2425]	1731
	1732
[3281]	1733	NAMELIST /chemistry_parameters/ bc_cs_b, &
[3085]	1734	bc_cs_t, &
	1735	call_chem_at_all_substeps, &
	1736	chem_debug0, &
	1737	chem_debug1, &
	1738	chem_debug2, &
	1739	chem_gasphase_on, &
	1740	cs_heights, &
	1741	cs_name, &
	1742	cs_profile, &
	1743	cs_profile_name, &
	1744	cs_surface, &
	1745	decycle_chem_lr, &
	1746	decycle_chem_ns, &
	1747	decycle_method, &
	1748	emiss_factor_main, &
	1749	emiss_factor_side, &
	1750	icntrl, &
	1751	main_street_id, &
	1752	max_street_id, &
	1753	my_steps, &
[3228]	1754	nest_chemistry, &
[3085]	1755	rcntrl, &
	1756	side_street_id, &
	1757	photolysis_scheme, &
	1758	wall_csflux, &
	1759	cs_vertical_gradient, &
	1760	top_csflux, &
	1761	surface_csflux, &
	1762	surface_csflux_name, &
	1763	cs_surface_initial_change, &
[3190]	1764	cs_vertical_gradient_level, &
[3282]	1765	! namelist parameters for emissions
[3190]	1766	mode_emis, &
	1767	time_fac_type, &
	1768	daytype_mdh, &
	1769	do_emis
[3085]	1770
	1771	!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
	1772	!-- so this way we could prescribe a specific flux value for each species
	1773	!> chemistry_parameters for initial profiles
	1774	!> cs_names = 'O3', 'NO2', 'NO', ... to set initial profiles)
	1775	!> cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
	1776	!> cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
	1777	!> cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, ..... (chem spcs conc at height lvls chem_heights(1,:)) etc.
	1778	!> If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
	1779	!> then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
[2425]	1780
	1781	!
[3085]	1782	!-- Read chem namelist
[3185]	1783
[3281]	1784	INTEGER :: ier
	1785	CHARACTER(LEN=64) :: text
	1786	CHARACTER(LEN=8) :: solver_type
[3185]	1787
[3281]	1788	icntrl = 0
	1789	rcntrl = 0.0_wp
	1790	my_steps = 0.0_wp
	1791	photolysis_scheme = 'simple'
	1792	atol = 1.0_wp
	1793	rtol = 0.01_wp
	1794	!
	1795	!-- Try to find chemistry package
	1796	REWIND ( 11 )
	1797	line = ' '
	1798	DO WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
[3287]	1799	READ ( 11, '(A)', END=20 ) line
[3281]	1800	ENDDO
	1801	BACKSPACE ( 11 )
	1802	!
	1803	!-- Read chemistry namelist
[3287]	1804	READ ( 11, chemistry_parameters, ERR = 10, END = 20 )
[3281]	1805	!
	1806	!-- Enable chemistry model
	1807	air_chemistry = .TRUE.
[3287]	1808	GOTO 20
	1809
	1810	10 BACKSPACE( 11 )
	1811	READ( 11 , '(A)') line
	1812	CALL parin_fail_message( 'chemistry_parameters', line )
	1813
	1814	20 CONTINUE
	1815
	1816	!
	1817	!-- check for emission mode for chem species
[3281]	1818	IF ( (mode_emis /= 'PARAMETERIZED') .AND. ( mode_emis /= 'DEFAULT') .AND. (mode_emis /= 'PRE-PROCESSED' ) ) THEN
	1819	message_string = 'Incorrect mode_emiss option select. Please check spelling'
	1820	CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
	1821	ENDIF
[2425]	1822
[3282]	1823	t_steps = my_steps
	1824
	1825	!-- Determine the number of user-defined profiles and append them to the
	1826	!-- standard data output (data_output_pr)
[3085]	1827	max_pr_cs_tmp = 0
	1828	i = 1
[2425]	1829
[3085]	1830	DO WHILE ( data_output_pr(i) /= ' ' .AND. i <= 100 )
	1831	IF ( TRIM( data_output_pr(i)(1:3)) == 'kc_' ) THEN
	1832	max_pr_cs_tmp = max_pr_cs_tmp+1
[2615]	1833	ENDIF
[3085]	1834	i = i +1
[2615]	1835	ENDDO
[2425]	1836
[3085]	1837	IF ( max_pr_cs_tmp > 0 ) THEN
	1838	cs_pr_namelist_found = .TRUE.
	1839	max_pr_cs = max_pr_cs_tmp
[2914]	1840	ENDIF
[3185]	1841
	1842	! Set Solver Type
[3260]	1843	IF(icntrl(3) == 0) THEN
[3185]	1844	solver_type = 'rodas3' !Default
[3260]	1845	ELSE IF(icntrl(3) == 1) THEN
[3185]	1846	solver_type = 'ros2'
[3260]	1847	ELSE IF(icntrl(3) == 2) THEN
[3185]	1848	solver_type = 'ros3'
[3260]	1849	ELSE IF(icntrl(3) == 3) THEN
[3185]	1850	solver_type = 'ro4'
[3260]	1851	ELSE IF(icntrl(3) == 4) THEN
[3185]	1852	solver_type = 'rodas3'
[3260]	1853	ELSE IF(icntrl(3) == 5) THEN
[3185]	1854	solver_type = 'rodas4'
[3260]	1855	ELSE IF(icntrl(3) == 6) THEN
[3185]	1856	solver_type = 'Rang3'
[3260]	1857	ELSE
[3373]	1858	message_string = 'illegal solver type'
	1859	CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
[3260]	1860	END IF
[3185]	1861
[3373]	1862	!
	1863	!-- todo: remove or replace by "CALL message" mechanism (kanani)
	1864	! write(text,*) 'gas_phase chemistry: solver_type = ',trim(solver_type)
[3185]	1865	!kk Has to be changed to right calling sequence
	1866	!kk CALL location_message( trim(text), .FALSE. )
[3373]	1867	! IF(myid == 0) THEN
	1868	! write(9,*) ' '
	1869	! write(9,*) 'kpp setup '
	1870	! write(9,*) ' '
	1871	! write(9,*) ' gas_phase chemistry: solver_type = ',trim(solver_type)
	1872	! write(9,*) ' '
	1873	! write(9,*) ' Hstart = ',rcntrl(3)
	1874	! write(9,*) ' FacMin = ',rcntrl(4)
	1875	! write(9,*) ' FacMax = ',rcntrl(5)
	1876	! write(9,*) ' '
	1877	! IF(vl_dim > 1) THEN
	1878	! write(9,*) ' Vector mode vektor length = ',vl_dim
	1879	! ELSE
	1880	! write(9,*) ' Scalar mode'
	1881	! ENDIF
	1882	! write(9,*) ' '
	1883	! END IF
[3185]	1884
[3281]	1885	RETURN
[2491]	1886
[3281]	1887	END SUBROUTINE chem_parin
[2467]	1888
	1889	!
[2828]	1890	!------------------------------------------------------------------------------!
[2467]	1891	!
[2828]	1892	! Description:
	1893	! ------------
	1894	!> Subroutine calculating prognostic equations for chemical species
	1895	!> (vector-optimized).
	1896	!> Routine is called separately for each chemical species over a loop from
	1897	!> prognostic_equations.
	1898	!------------------------------------------------------------------------------!
[3281]	1899	SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m, &
	1900	pr_init_cs, ilsp )
[2467]	1901
[3281]	1902	USE advec_s_pw_mod, &
	1903	ONLY: advec_s_pw
	1904	USE advec_s_up_mod, &
	1905	ONLY: advec_s_up
	1906	USE advec_ws, &
	1907	ONLY: advec_s_ws
	1908	USE diffusion_s_mod, &
	1909	ONLY: diffusion_s
	1910	USE indices, &
	1911	ONLY: nxl, nxr, nyn, nys, wall_flags_0
	1912	USE pegrid
	1913	USE surface_mod, &
	1914	ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	1915	surf_usm_v
[2828]	1916
[3281]	1917	IMPLICIT NONE
[2828]	1918
[3281]	1919	INTEGER :: i !< running index
	1920	INTEGER :: j !< running index
	1921	INTEGER :: k !< running index
[2828]	1922
[3281]	1923	INTEGER(iwp),INTENT(IN) :: ilsp !<
[2828]	1924
[3281]	1925	REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !<
[2828]	1926
[3281]	1927	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar !<
	1928	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p !<
	1929	REAL(wp), DIMENSION(:,:,:), POINTER :: tcs_scalar_m !<
[2828]	1930
	1931
[3281]	1932	!
	1933	!-- Tendency terms for chemical species
	1934	tend = 0.0_wp
	1935	!
	1936	!-- Advection terms
	1937	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1938	IF ( ws_scheme_sca ) THEN
	1939	CALL advec_s_ws( cs_scalar, 'kc' )
	1940	ELSE
	1941	CALL advec_s_pw( cs_scalar )
	1942	ENDIF
[2828]	1943	ELSE
[3281]	1944	CALL advec_s_up( cs_scalar )
[2828]	1945	ENDIF
[3281]	1946	!
	1947	!-- Diffusion terms (the last three arguments are zero)
	1948	CALL diffusion_s( cs_scalar, &
	1949	surf_def_h(0)%cssws(ilsp,:), &
	1950	surf_def_h(1)%cssws(ilsp,:), &
	1951	surf_def_h(2)%cssws(ilsp,:), &
	1952	surf_lsm_h%cssws(ilsp,:), &
	1953	surf_usm_h%cssws(ilsp,:), &
	1954	surf_def_v(0)%cssws(ilsp,:), &
	1955	surf_def_v(1)%cssws(ilsp,:), &
	1956	surf_def_v(2)%cssws(ilsp,:), &
	1957	surf_def_v(3)%cssws(ilsp,:), &
	1958	surf_lsm_v(0)%cssws(ilsp,:), &
	1959	surf_lsm_v(1)%cssws(ilsp,:), &
	1960	surf_lsm_v(2)%cssws(ilsp,:), &
	1961	surf_lsm_v(3)%cssws(ilsp,:), &
	1962	surf_usm_v(0)%cssws(ilsp,:), &
	1963	surf_usm_v(1)%cssws(ilsp,:), &
	1964	surf_usm_v(2)%cssws(ilsp,:), &
	1965	surf_usm_v(3)%cssws(ilsp,:) )
	1966	!
	1967	!-- Prognostic equation for chemical species
	1968	DO i = nxl, nxr
	1969	DO j = nys, nyn
	1970	DO k = nzb+1, nzt
	1971	cs_scalar_p(k,j,i) = cs_scalar(k,j,i) &
	1972	+ ( dt_3d * &
	1973	( tsc(2) * tend(k,j,i) &
	1974	+ tsc(3) * tcs_scalar_m(k,j,i) &
	1975	) &
	1976	- tsc(5) * rdf_sc(k) &
	1977	* ( cs_scalar(k,j,i) - pr_init_cs(k) ) &
	1978	) &
	1979	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2828]	1980
[3281]	1981	IF ( cs_scalar_p(k,j,i) < 0.0_wp ) cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)
	1982	ENDDO
[2828]	1983	ENDDO
	1984	ENDDO
[3281]	1985	!
	1986	!-- Calculate tendencies for the next Runge-Kutta step
	1987	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1988	IF ( intermediate_timestep_count == 1 ) THEN
	1989	DO i = nxl, nxr
	1990	DO j = nys, nyn
	1991	DO k = nzb+1, nzt
	1992	tcs_scalar_m(k,j,i) = tend(k,j,i)
	1993	ENDDO
[2828]	1994	ENDDO
	1995	ENDDO
[3281]	1996	ELSEIF ( intermediate_timestep_count < &
	1997	intermediate_timestep_count_max ) THEN
	1998	DO i = nxl, nxr
	1999	DO j = nys, nyn
	2000	DO k = nzb+1, nzt
	2001	tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i) &
	2002	+ 5.3125_wp * tcs_scalar_m(k,j,i)
	2003	ENDDO
[2828]	2004	ENDDO
	2005	ENDDO
[3281]	2006	ENDIF
[2828]	2007	ENDIF
	2008
[3281]	2009	END SUBROUTINE chem_prognostic_equations
[3228]	2010
[2482]	2011	!------------------------------------------------------------------------------!
[2535]	2012	!
[2482]	2013	! Description:
	2014	! ------------
[3085]	2015	!> Subroutine calculating prognostic equations for chemical species
	2016	!> (cache-optimized).
	2017	!> Routine is called separately for each chemical species over a loop from
	2018	!> prognostic_equations.
[2482]	2019	!------------------------------------------------------------------------------!
[3281]	2020	SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs, &
	2021	i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, &
	2022	flux_l_cs, diss_l_cs )
	2023	USE pegrid
	2024	USE advec_ws, ONLY: advec_s_ws
	2025	USE advec_s_pw_mod, ONLY: advec_s_pw
	2026	USE advec_s_up_mod, ONLY: advec_s_up
	2027	USE diffusion_s_mod, ONLY: diffusion_s
	2028	USE indices, ONLY: wall_flags_0
	2029	USE surface_mod, ONLY: surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	2030	surf_usm_v
[3085]	2031
	2032
[3281]	2033	IMPLICIT NONE
[3085]	2034
[3281]	2035	REAL(wp), DIMENSION(:,:,:), POINTER :: cs_scalar_p, cs_scalar, tcs_scalar_m
[3085]	2036
[3281]	2037	INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp
	2038	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: flux_s_cs !<
	2039	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: diss_s_cs !<
	2040	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: flux_l_cs !<
	2041	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: diss_l_cs !<
	2042	REAL(wp), DIMENSION(0:nz+1) :: pr_init_cs !<
[3085]	2043
	2044	!-- local variables
	2045
[3281]	2046	INTEGER :: k
	2047	!
	2048	!-- Tendency-terms for chem spcs.
	2049	tend(:,j,i) = 0.0_wp
[3085]	2050	!
	2051	!-- Advection terms
[3281]	2052	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2053	IF ( ws_scheme_sca ) THEN
	2054	CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs, &
	2055	flux_l_cs, diss_l_cs, i_omp_start, tn )
	2056	ELSE
	2057	CALL advec_s_pw( i, j, cs_scalar )
	2058	ENDIF
[3085]	2059	ELSE
[3281]	2060	CALL advec_s_up( i, j, cs_scalar )
[3085]	2061	ENDIF
	2062
	2063	!
	2064
	2065	!-- Diffusion terms (the last three arguments are zero)
	2066
[3281]	2067	CALL diffusion_s( i, j, cs_scalar, &
	2068	surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:), &
	2069	surf_def_h(2)%cssws(ilsp,:), &
	2070	surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:), &
	2071	surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:), &
	2072	surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:), &
	2073	surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:), &
	2074	surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:), &
	2075	surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:), &
	2076	surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
	2077
[3085]	2078	!
	2079	!-- Prognostic equation for chem spcs
[3281]	2080	DO k = nzb+1, nzt
	2081	cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d * &
	2082	( tsc(2) * tend(k,j,i) + &
	2083	tsc(3) * tcs_scalar_m(k,j,i) ) &
	2084	- tsc(5) * rdf_sc(k) &
	2085	* ( cs_scalar(k,j,i) - pr_init_cs(k) ) &
	2086	) &
	2087	* MERGE( 1.0_wp, 0.0_wp, &
	2088	BTEST( wall_flags_0(k,j,i), 0 ) &
	2089	)
[2467]	2090
[3281]	2091	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
	2092	ENDDO
[3085]	2093
[2482]	2094	!
[3085]	2095	!-- Calculate tendencies for the next Runge-Kutta step
[3281]	2096	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2097	IF ( intermediate_timestep_count == 1 ) THEN
	2098	DO k = nzb+1, nzt
	2099	tcs_scalar_m(k,j,i) = tend(k,j,i)
	2100	ENDDO
	2101	ELSEIF ( intermediate_timestep_count < &
	2102	intermediate_timestep_count_max ) THEN
	2103	DO k = nzb+1, nzt
	2104	tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	2105	5.3125_wp * tcs_scalar_m(k,j,i)
	2106	ENDDO
	2107	ENDIF
[3085]	2108	ENDIF
	2109
[3281]	2110	END SUBROUTINE chem_prognostic_equations_ij
[3085]	2111
[2482]	2112	!
[3085]	2113	!------------------------------------------------------------------------------!
	2114	!
	2115	! Description:
	2116	! ------------
	2117	!> Subroutine to read restart data of chemical species
	2118	!------------------------------------------------------------------------------!
[2828]	2119
[3281]	2120	SUBROUTINE chem_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
	2121	nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc, &
	2122	nys_on_file, tmp_3d, found )
	2123
	2124	USE control_parameters
	2125
	2126	USE indices
	2127
	2128	USE pegrid
[3085]	2129
[3281]	2130	IMPLICIT NONE
[3085]	2131
[3281]	2132	CHARACTER (LEN=20) :: spc_name_av !<
	2133
	2134	INTEGER(iwp) :: i, lsp !<
	2135	INTEGER(iwp) :: k !<
	2136	INTEGER(iwp) :: nxlc !<
	2137	INTEGER(iwp) :: nxlf !<
	2138	INTEGER(iwp) :: nxl_on_file !<
	2139	INTEGER(iwp) :: nxrc !<
	2140	INTEGER(iwp) :: nxrf !<
	2141	INTEGER(iwp) :: nxr_on_file !<
	2142	INTEGER(iwp) :: nync !<
	2143	INTEGER(iwp) :: nynf !<
	2144	INTEGER(iwp) :: nyn_on_file !<
	2145	INTEGER(iwp) :: nysc !<
	2146	INTEGER(iwp) :: nysf !<
	2147	INTEGER(iwp) :: nys_on_file !<
[3085]	2148
[3281]	2149	LOGICAL, INTENT(OUT) :: found
[3085]	2150
[3281]	2151	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
[3085]	2152
	2153
[3281]	2154	found = .FALSE.
[3085]	2155
	2156
[3281]	2157	IF ( ALLOCATED(chem_species) ) THEN
[3085]	2158
[3281]	2159	DO lsp = 1, nspec
[3085]	2160
[3281]	2161	!< for time-averaged chemical conc.
	2162	spc_name_av = TRIM(chem_species(lsp)%name)//'_av'
[3085]	2163
[3281]	2164	IF (restart_string(1:length) == TRIM(chem_species(lsp)%name) ) &
	2165	THEN
	2166	!< read data into tmp_3d
	2167	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2168	!< fill ..%conc in the restart run
	2169	chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp, &
	2170	nxlc-nbgp:nxrc+nbgp) = &
	2171	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2172	found = .TRUE.
	2173	ELSEIF (restart_string(1:length) == spc_name_av ) THEN
	2174	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2175	chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp, &
	2176	nxlc-nbgp:nxrc+nbgp) = &
	2177	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2178	found = .TRUE.
	2179	ENDIF
[3085]	2180
[3281]	2181	ENDDO
[3085]	2182
[3281]	2183	ENDIF
[3085]	2184
	2185
[3281]	2186	END SUBROUTINE chem_rrd_local
[3085]	2187	!
	2188	!-------------------------------------------------------------------------------!
	2189	!> Description:
	2190	!> Calculation of horizontally averaged profiles
	2191	!> This routine is called for every statistic region (sr) defined by the user,
	2192	!> but at least for the region "total domain" (sr=0).
	2193	!> quantities.
[2828]	2194	!------------------------------------------------------------------------------!
[3281]	2195	SUBROUTINE chem_statistics( mode, sr, tn )
[3085]	2196
	2197
[3281]	2198	USE arrays_3d
	2199	USE indices
	2200	USE kinds
	2201	USE pegrid
	2202	USE statistics
[3085]	2203
	2204	USE user
	2205
	2206	IMPLICIT NONE
	2207
	2208	CHARACTER (LEN=*) :: mode !<
	2209
	2210	INTEGER(iwp) :: i !< running index on x-axis
	2211	INTEGER(iwp) :: j !< running index on y-axis
	2212	INTEGER(iwp) :: k !< vertical index counter
	2213	INTEGER(iwp) :: sr !< statistical region
	2214	INTEGER(iwp) :: tn !< thread number
	2215	INTEGER(iwp) :: n !<
	2216	INTEGER(iwp) :: m !<
	2217	INTEGER(iwp) :: lpr !< running index chem spcs
	2218	! REAL(wp), &
	2219	! DIMENSION(dots_num_palm+1:dots_max) :: &
	2220	! ts_value_l !<
	2221
	2222	IF ( mode == 'profiles' ) THEN
[2615]	2223	!
[3085]	2224	!-- Sample on how to calculate horizontally averaged profiles of user-
	2225	!-- defined quantities. Each quantity is identified by the index
	2226	!-- "pr_palm+#" where "#" is an integer starting from 1. These
	2227	!-- user-profile-numbers must also be assigned to the respective strings
	2228	!-- given by data_output_pr_cs in routine user_check_data_output_pr.
	2229	!-- hom(:,:,:,:) = dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
	2230	! wpt), dim-4 = statistical region.
	2231
	2232	!$OMP DO
	2233	DO i = nxl, nxr
	2234	DO j = nys, nyn
	2235	DO k = nzb, nzt+1
	2236	DO lpr = 1, cs_pr_count
	2237
	2238	sums_l(k,pr_palm+lpr,tn) = sums_l(k,pr_palm+lpr,tn) + &
	2239	chem_species(cs_pr_index(lpr))%conc(k,j,i) * &
	2240	rmask(j,i,sr) * &
	2241	MERGE( 1.0_wp, 0.0_wp, &
	2242	BTEST( wall_flags_0(k,j,i), 22 ) )
	2243	ENDDO
	2244	ENDDO
	2245	ENDDO
	2246	ENDDO
	2247	ELSEIF ( mode == 'time_series' ) THEN
[3282]	2248	CALL location_message( 'Time series not calculated for chemistry', .TRUE. )
	2249	ENDIF
[3085]	2250
[3281]	2251	END SUBROUTINE chem_statistics
[3085]	2252	!
	2253	!------------------------------------------------------------------------------!
	2254	!
[2828]	2255	! Description:
	2256	! ------------
[3085]	2257	!> Subroutine for swapping of timelevels for chemical species
	2258	!> called out from subroutine swap_timelevel
	2259	!------------------------------------------------------------------------------!
	2260
[3281]	2261	SUBROUTINE chem_swap_timelevel (level)
[3085]	2262
[3281]	2263	IMPLICIT none
[3085]	2264
[3281]	2265	INTEGER,INTENT(IN) :: level
	2266	!-- local variables
	2267	INTEGER :: lsp
[3085]	2268
[3298]	2269
[3281]	2270	IF ( level == 0 ) THEN
	2271	DO lsp=1, nvar
	2272	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp)
	2273	chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp)
	2274	ENDDO
	2275	ELSE
	2276	DO lsp=1, nvar
	2277	chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_2(:,:,:,lsp)
	2278	chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) => spec_conc_1(:,:,:,lsp)
	2279	ENDDO
	2280	ENDIF
[3085]	2281
[3281]	2282	RETURN
	2283	END SUBROUTINE chem_swap_timelevel
[3085]	2284	!
[2615]	2285	!------------------------------------------------------------------------------!
[3085]	2286	!
[2615]	2287	! Description:
	2288	! ------------
[3085]	2289	!> Subroutine to write restart data for chemistry model
[2615]	2290	!------------------------------------------------------------------------------!
[3281]	2291	SUBROUTINE chem_wrd_local
[2482]	2292
[3281]	2293	IMPLICIT NONE
[2615]	2294
[3281]	2295	INTEGER(iwp) :: lsp !<
[2615]	2296
[3281]	2297	DO lsp = 1, nspec
	2298	CALL wrd_write_string( TRIM( chem_species(lsp)%name ))
	2299	WRITE ( 14 ) chem_species(lsp)%conc
	2300	CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
	2301	WRITE ( 14 ) chem_species(lsp)%conc_av
	2302	ENDDO
[2615]	2303
[3281]	2304	END SUBROUTINE chem_wrd_local
[2615]	2305	END MODULE chemistry_model_mod
	2306

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |