Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |