Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |