Home

Context Navigation

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

Last change on this file since 2951 was 2932, checked in by maronga, 7 years ago
renamed all Fortran NAMELISTS
Property svn:keywords set to `Id`
File size: 74.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |