Home

Context Navigation

source: palm/trunk/SOURCE/chem_emissions_mod.f90 @ 4661

Last change on this file since 4661 was 4559, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 97.6 KB

Line
1	!> @file chem_emissions_mod.f90
2	!--------------------------------------------------------------------------------------------------!
3	! This file is part of PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
7	! (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11	! Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
14	! <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 2018-2020 Leibniz Universitaet Hannover
17	! Copyright 2018-2020 Freie Universitaet Berlin
18	! Copyright 2018-2020 Karlsruhe Institute of Technology
19	!--------------------------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! ------------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: chem_emissions_mod.f90 4559 2020-06-11 08:51:48Z pavelkrc $
28	! file re-formatted to follow the PALM coding standard
29	!
30	! 4481 2020-03-31 18:55:54Z maronga
31	! Implemented on-demand read mode for LOD 2 (ECC)
32	! - added following module global variables
33	! - input_file_chem (namesake in netcdf_data_input_mod is local)
34	! - timestamps
35	! - added following public subroutines / interfaces
36	! - chem_emisisons_header_init
37	! - chem_emisisons_update_on_demand
38	! - added following local subroutines
39	! - chem_emisisons_header_init_lod2
40	! - chem_emisisons_update_on_demand_lod2
41	! - added following local auxiliary subroutines
42	! - chem_emissions_init_species ( )
43	! - chem_emissions_init_timestamps ( )
44	! - chem_emissions_assign_surface_flux ( )
45	! - added following local functions
46	! - chem_emisisons_convert_base_units ( )
47	! - chem_emissions_mass_2_molar_flux ( )
48	! - chem_emissions_locate_species ( )
49	! - chem_emissions_locate_timestep ( )
50	! - added following error messages
51	! - CM0468 - LOD mismatch (namelist / chemistry file)
52	! - CM0469 - Timestamps no in choronological order
53	! - depreciated unused module variable filename_emis
54	!
55	! 4356 2019-12-20 17:09:33Z suehring
56	! Minor formatting adjustment
57	!
58	! 4242 2019-09-27 12:59:10Z suehring
59	! Adjust index_hh access to new definition accompanied with new
60	! palm_date_time_mod. Note, this is just a preleminary fix. (E Chan)
61	!
62	! 4230 2019-09-11 13:58:14Z suehring
63	! Bugfix, consider that time_since_reference_point can be also negative when
64	! time indices are determined.
65	!
66	! 4227 2019-09-10 18:04:34Z gronemeier
67	! implement new palm_date_time_mod
68	!
69	! 4223 2019-09-10 09:20:47Z gronemeier
70	! Unused routine chem_emissions_check_parameters commented out due to uninitialized content
71	!
72	! 4182 2019-08-22 15:20:23Z scharf
73	! Corrected "Former revisions" section
74	!
75	! 4157 2019-08-14 09:19:12Z suehring
76	! Replace global arrays also in mode_emis branch
77	!
78	! 4154 2019-08-13 13:35:59Z suehring
79	! Replace global arrays for emissions by local ones.
80	!
81	! 4144 2019-08-06 09:11:47Z raasch
82	! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
83	!
84	! 4055 2019-06-27 09:47:29Z suehring
85	! - replaced local thermo. constants w/ module definitions in basic_constants_and_equations_mod
86	! (rgas_univ, p_0, r_d_cp)
87	! - initialize array emis_distribution immediately following allocation
88	! - lots of minor formatting changes based on review sesson in 20190325 (E.C. Chan)
89	!
90	! 3968 2019-05-13 11:04:01Z suehring
91	! - in subroutine chem_emissions_match replace all decision structures relating to mode_emis to
92	! emiss_lod
93	! - in subroutine chem_check_parameters replace emt%nspec with emt%n_emiss_species
94	! - spring cleaning (E.C. Chan)
95	!
96	! 3885 2019-04-11 11:29:34Z kanani
97	! Changes related to global restructuring of location messages and introduction of additional debug
98	! messages
99	!
100	! 3831 2019-03-28 09:11:22Z forkel
101	! added nvar to USE chem_gasphase_mod (chem_modules will not include nvar anymore)
102	!
103	! 3788 2019-03-07 11:40:09Z banzhafs
104	! Removed unused variables from chem_emissions_mod
105	!
106	! 3772 2019-02-28 15:51:57Z suehring
107	! - In case of parametrized emissions, assure that emissions are only on natural surfaces
108	! (i.e. streets) and not on urban surfaces.
109	! - some unnecessary if clauses removed
110	!
111	! 3685 2019 -01-21 01:02:11Z knoop
112	! Some interface calls moved to module_interface + cleanup
113	! 3286 2018-09-28 07:41:39Z forkel
114	!
115	! Authors:
116	! --------
117	! @author Emmanuele Russo (FU-Berlin)
118	! @author Sabine Banzhaf (FU-Berlin)
119	! @author Martijn Schaap (FU-Berlin, TNO Utrecht)
120	!
121	! Description:
122	! ------------
123	!> MODULE for reading-in Chemistry Emissions data
124	!>
125	!> @todo Rename nspec to n_emis to avoid inteferece with nspec from chem_gasphase_mod
126	!> @todo Check_parameters routine should be developed: for now it includes just one condition
127	!> @todo Use of Restart files not contempled at the moment
128	!> @todo revise indices of files read from the netcdf: preproc_emission_data and expert_emission_data
129	!> @todo for now emission data may be passed on a singular vertical level: need to be more flexible
130	!> @todo fill/activate restart option in chem_emissions_init
131	!> @todo discuss dt_emis
132	!> @note <Enter notes on the module>
133	!> @bug <Enter known bugs here>
134	!--------------------------------------------------------------------------------------------------!
135
136	MODULE chem_emissions_mod
137
138	USE arrays_3d, &
139	ONLY: rho_air
140
141	USE basic_constants_and_equations_mod, &
142	ONLY: p_0, rd_d_cp, rgas_univ
143
144	USE control_parameters, &
145	ONLY: debug_output, end_time, initializing_actions, intermediate_timestep_count, &
146	message_string, dt_3d
147
148	USE indices
149
150	USE kinds
151
152	#if defined( __netcdf )
153	USE netcdf
154	#endif
155
156	USE netcdf_data_input_mod, &
157	ONLY: chem_emis_att_type, chem_emis_val_type
158
159	USE chem_gasphase_mod, &
160	ONLY: nvar, spc_names
161
162	USE chem_modules
163
164	USE statistics, &
165	ONLY: weight_pres
166
167	!
168	!-- 20200203 (ECC)
169	!-- Added new palm_date_time_mod for on-demand emission reading
170
171	USE palm_date_time_mod, &
172	ONLY: get_date_time
173
174	IMPLICIT NONE
175
176	!
177	!-- Declare all global variables within the module
178
179	!
180	!-- 20200203 (ECC) - variable unused
181	! CHARACTER (LEN=80) :: filename_emis !< Variable for the name of the netcdf input file
182
183	!
184	!-- 20200203 (ECC) new variables for on-demand read mode
185
186	CHARACTER(LEN=*), PARAMETER :: input_file_chem = 'PIDS_CHEM' !< chemistry file
187
188	CHARACTER(LEN=512), ALLOCATABLE, DIMENSION(:) :: timestamps !< timestamps in chemistry file
189
190
191	INTEGER(iwp) :: dt_emis !< Time Step Emissions
192	INTEGER(iwp) :: i !< index 1st selected dimension (some dims are not spatial)
193	INTEGER(iwp) :: j !< index 2nd selected dimension
194	INTEGER(iwp) :: i_end !< Index to end read variable from netcdf in one dims
195	INTEGER(iwp) :: i_start !< Index to start read variable from netcdf along one dims
196	INTEGER(iwp) :: j_end !< Index to end read variable from netcdf in additional dims
197	INTEGER(iwp) :: j_start !< Index to start read variable from netcdf in additional dims
198	INTEGER(iwp) :: len_index !< length of index (used for several indices)
199	INTEGER(iwp) :: len_index_pm !< length of PMs index
200	INTEGER(iwp) :: len_index_voc !< length of voc index
201	INTEGER(iwp) :: previous_timestamp_index !< index for current timestamp (20200203 ECC)
202	INTEGER(iwp) :: z_end !< Index to end read variable from netcdf in additional dims
203	INTEGER(iwp) :: z_start !< Index to start read variable from netcdf in additional dims
204
205	REAL(wp) :: conversion_factor !< Units Conversion Factor
206
207	SAVE
208
209	! !
210	! !-- Checks Input parameters
211	! INTERFACE chem_emissions_check_parameters
212	! MODULE PROCEDURE chem_emissions_check_parameters
213	! END INTERFACE chem_emissions_check_parameters
214	!
215	!-- Matching Emissions actions
216	INTERFACE chem_emissions_match
217	MODULE PROCEDURE chem_emissions_match
218	END INTERFACE chem_emissions_match
219	!
220	!-- Initialization actions
221	INTERFACE chem_emissions_init
222	MODULE PROCEDURE chem_emissions_init
223	END INTERFACE chem_emissions_init
224	!
225	!-- Setup of Emissions
226	INTERFACE chem_emissions_setup
227	MODULE PROCEDURE chem_emissions_setup
228	END INTERFACE chem_emissions_setup
229
230	!
231	!-- 20200203 (ECC) new interfaces for on-demand mode
232
233	!
234	!-- initialization actions for on-demand mode
235	INTERFACE chem_emissions_header_init
236	MODULE PROCEDURE chem_emissions_header_init
237	END INTERFACE chem_emissions_header_init
238	!
239	!-- load emission data for on-demand mode
240	INTERFACE chem_emissions_update_on_demand
241	MODULE PROCEDURE chem_emissions_update_on_demand
242	END INTERFACE chem_emissions_update_on_demand
243
244	!
245	!-- 20200203 (ECC) update public routines
246
247	! PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup
248
249	PUBLIC chem_emissions_init, chem_emissions_match, chem_emissions_setup, &
250	chem_emissions_header_init, chem_emissions_update_on_demand
251	!
252	!-- Public Variables
253	PUBLIC conversion_factor, len_index, len_index_pm, len_index_voc
254
255	CONTAINS
256
257	! !------------------------------------------------------------------------------------------------!
258	! ! Description:
259	! ! ------------
260	! !> Routine for checking input parameters
261	! !------------------------------------------------------------------------------------------------!
262	! SUBROUTINE chem_emissions_check_parameters
263	!
264	! IMPLICIT NONE
265	!
266	! TYPE(chem_emis_att_type) :: emt
267	!
268	! !
269	! !-- Check if species count matches the number of names
270	! !-- passed for the chemiscal species
271	!
272	! IF ( SIZE(emt%species_name) /= emt%n_emiss_species ) THEN
273	! ! IF ( SIZE(emt%species_name) /= emt%nspec ) THEN
274	!
275	! message_string = 'Numbers of input emission species names and number of species' // &
276	! 'for which emission values are given do not match'
277	! CALL message( 'chem_emissions_check_parameters', 'CM0437', 2, 2, 0, 6, 0 )
278	!
279	! ENDIF
280	!
281	! END SUBROUTINE chem_emissions_check_parameters
282
283
284	!--------------------------------------------------------------------------------------------------!
285	! Description:
286	! ------------
287	!> Matching the chemical species indices. The routine checks what are the indices of the emission
288	!> input species and the corresponding ones of the model species. The routine gives as output a
289	!> vector containing the number of common species: it is important to note that while the model
290	!> species are distinct, their values could be given to a single species in input.
291	!> For example, in the case of NO2 and NO, values may be passed in input as NOX values.
292	!--------------------------------------------------------------------------------------------------!
293
294	SUBROUTINE chem_emissions_match( emt_att,len_index )
295
296	INTEGER(iwp) :: ind_inp !< Parameters for cycling through chemical input species
297	INTEGER(iwp) :: ind_mod !< Parameters for cycling through chemical model species
298	INTEGER(iwp) :: ind_voc !< Indices to check whether a split for voc should be done
299	INTEGER(iwp) :: ispec !< index for cycle over effective number of emission species
300	INTEGER(iwp) :: nspec_emis_inp !< Variable where to store # of emission species in input
301
302	INTEGER(iwp), INTENT(INOUT) :: len_index !< number of common species between input dataset & model species
303
304	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< Chemistry Emission Array (decl. netcdf_data_input.f90)
305
306
307	IF ( debug_output ) CALL debug_message( 'chem_emissions_match', 'start' )
308
309	!
310	!-- Number of input emission species
311
312	nspec_emis_inp = emt_att%n_emiss_species
313	! nspec_emis_inp=emt_att%nspec
314
315	!
316	!-- Check the emission LOD: 0 (PARAMETERIZED), 1 (DEFAULT), 2 (PRE-PROCESSED)
317	SELECT CASE (emiss_lod)
318
319	!
320	!-- LOD 0 (PARAMETERIZED mode)
321	CASE (0)
322
323	len_index = 0
324	!
325	!-- Number of species and number of matched species can be different but call is only made if
326	!-- both are greater than zero.
327	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
328
329	!
330	!-- Cycle over model species
331	DO ind_mod = 1, nvar
332	ind_inp = 1
333	DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' ) !< 'novalue' is the default
334	IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
335	len_index = len_index + 1
336	ENDIF
337	ind_inp = ind_inp + 1
338	ENDDO
339	ENDDO
340
341	IF ( len_index > 0 ) THEN
342
343	!
344	!-- Allocation of Arrays of the matched species
345	ALLOCATE ( match_spec_input(len_index) )
346	ALLOCATE ( match_spec_model(len_index) )
347
348	!
349	!-- Pass species indices to declared arrays
350	len_index = 0
351
352	DO ind_mod = 1, nvar
353	ind_inp = 1
354	DO WHILE ( TRIM( surface_csflux_name(ind_inp) ) /= 'novalue' )
355	IF ( TRIM( surface_csflux_name(ind_inp) ) == TRIM(spc_names(ind_mod) ) ) &
356	THEN
357	len_index = len_index + 1
358	match_spec_input(len_index) = ind_inp
359	match_spec_model(len_index) = ind_mod
360	ENDIF
361	ind_inp = ind_inp + 1
362	END DO
363	END DO
364
365	!
366	!-- Check
367	DO ispec = 1, len_index
368
369	IF ( emiss_factor_main(match_spec_input(ispec) ) < 0 .AND. &
370	emiss_factor_side(match_spec_input(ispec) ) < 0 ) THEN
371
372	message_string = 'PARAMETERIZED emissions mode selected:' // &
373	' EMISSIONS POSSIBLE ONLY ON STREET SURFACES' // &
374	' but values of scaling factors for street types' // &
375	' emiss_factor_main AND emiss_factor_side' // &
376	' not provided for each of the emissions species' // &
377	' or not provided at all: PLEASE set a finite value' // &
378	' for these parameters in the chemistry namelist'
379	CALL message( 'chem_emissions_matching', 'CM0442', 2, 2, 0, 6, 0 )
380
381	ENDIF
382
383	END DO
384
385
386	ELSE
387
388	message_string = 'Non of given Emission Species' // &
389	' matches' // &
390	' model chemical species' // &
391	' Emission routine is not called'
392	CALL message( 'chem_emissions_matching', 'CM0443', 0, 0, 0, 6, 0 )
393
394	ENDIF
395
396	ELSE
397
398	message_string = 'Array of Emission species not allocated: ' // &
399	' Either no emission species are provided as input or' // &
400	' no chemical species are used by PALM.' // &
401	' Emission routine is not called'
402	CALL message( 'chem_emissions_matching', 'CM0444', 0, 2, 0, 6, 0 )
403
404	ENDIF
405
406	!
407	!-- LOD 1 (DEFAULT mode)
408	CASE (1)
409
410	len_index = 0 ! total number of species (to be accumulated)
411	len_index_voc = 0 ! total number of VOCs (to be accumulated)
412	len_index_pm = 3 ! total number of PMs: PM1, PM2.5, PM10.
413
414	!
415	!-- Number of model species and input species could be different but process this only when both are
416	!-- non-zero
417	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
418
419	!
420	!-- Cycle over model species
421	DO ind_mod = 1, nvar
422
423	!
424	!-- Cycle over input species
425	DO ind_inp = 1, nspec_emis_inp
426
427	!
428	!-- Check for VOC Species
429	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" ) THEN
430	DO ind_voc= 1, emt_att%nvoc
431
432	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) &
433	THEN
434	len_index = len_index + 1
435	len_index_voc = len_index_voc + 1
436	ENDIF
437
438	END DO
439	ENDIF
440
441	!
442	!-- PMs: There is one input species name for all PM. This variable has 3 dimensions, one for PM1,
443	!-- PM2.5 and PM10
444	IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" ) THEN
445
446	IF ( TRIM( spc_names(ind_mod) ) == "PM1" ) THEN
447	len_index = len_index + 1
448	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" ) THEN
449	len_index = len_index + 1
450	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" ) THEN
451	len_index = len_index + 1
452	ENDIF
453
454	ENDIF
455
456	!
457	!-- NOX: NO2 and NO
458	IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" ) THEN
459
460	IF ( TRIM( spc_names(ind_mod) ) == "NO" ) THEN
461	len_index = len_index + 1
462	ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" ) THEN
463	len_index = len_index + 1
464	ENDIF
465
466	ENDIF
467
468	!
469	!-- SOX: SO2 and SO4
470	IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
471
472	IF ( TRIM( spc_names(ind_mod) ) == "SO2" ) THEN
473	len_index = len_index + 1
474	ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" ) THEN
475	len_index = len_index + 1
476	ENDIF
477
478	ENDIF
479
480	!
481	!-- Other Species
482	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
483	len_index = len_index + 1
484	ENDIF
485
486	END DO ! ind_inp ...
487
488	END DO ! ind_mod ...
489
490
491	!
492	!-- Allocate arrays
493	IF ( len_index > 0 ) THEN
494
495	ALLOCATE ( match_spec_input(len_index) )
496	ALLOCATE ( match_spec_model(len_index) )
497
498	IF ( len_index_voc > 0 ) THEN
499	!
500	!-- Contains indices of the VOC model species
501	ALLOCATE( match_spec_voc_model(len_index_voc) )
502	!
503	!-- Contains the indices of different values of VOC composition of input variable
504	!-- VOC_composition
505	ALLOCATE( match_spec_voc_input(len_index_voc) )
506
507	ENDIF
508
509	!
510	!-- Pass the species indices to declared arrays
511	len_index = 0
512	len_index_voc = 0
513
514	DO ind_mod = 1, nvar
515	DO ind_inp = 1, nspec_emis_inp
516
517	!
518	!-- VOCs
519	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND. &
520	ALLOCATED( match_spec_voc_input ) ) THEN
521
522	DO ind_voc = 1, emt_att%nvoc
523
524	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )&
525	THEN
526
527	len_index = len_index + 1
528	len_index_voc = len_index_voc + 1
529
530	match_spec_input(len_index) = ind_inp
531	match_spec_model(len_index) = ind_mod
532
533	match_spec_voc_input(len_index_voc) = ind_voc
534	match_spec_voc_model(len_index_voc) = ind_mod
535
536	ENDIF
537
538	END DO
539
540	ENDIF
541
542	!
543	!-- PMs
544	IF ( TRIM( emt_att%species_name(ind_inp) ) == "PM" ) THEN
545
546	IF ( TRIM( spc_names(ind_mod) ) == "PM1" ) THEN
547	len_index = len_index + 1
548	match_spec_input(len_index) = ind_inp
549	match_spec_model(len_index) = ind_mod
550	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM25" ) THEN
551	len_index = len_index + 1
552	match_spec_input(len_index) = ind_inp
553	match_spec_model(len_index) = ind_mod
554	ELSEIF ( TRIM( spc_names(ind_mod) ) == "PM10" ) THEN
555	len_index = len_index + 1
556	match_spec_input(len_index) = ind_inp
557	match_spec_model(len_index) = ind_mod
558	ENDIF
559
560	ENDIF
561
562	!
563	!-- NOX
564	IF ( TRIM( emt_att%species_name(ind_inp) ) == "NOX" ) THEN
565
566	IF ( TRIM( spc_names(ind_mod) ) == "NO" ) THEN
567	len_index = len_index + 1
568
569	match_spec_input(len_index) = ind_inp
570	match_spec_model(len_index) = ind_mod
571
572	ELSEIF ( TRIM( spc_names(ind_mod) ) == "NO2" ) THEN
573	len_index = len_index + 1
574
575	match_spec_input(len_index) = ind_inp
576	match_spec_model(len_index) = ind_mod
577
578	ENDIF
579
580	ENDIF
581
582
583	!
584	!-- SOX
585	IF ( TRIM( emt_att%species_name(ind_inp) ) == "SOX" ) THEN
586
587	IF ( TRIM( spc_names(ind_mod) ) == "SO2" ) THEN
588	len_index = len_index + 1
589	match_spec_input(len_index) = ind_inp
590	match_spec_model(len_index) = ind_mod
591	ELSEIF ( TRIM( spc_names(ind_mod) ) == "SO4" ) THEN
592	len_index = len_index + 1
593	match_spec_input(len_index) = ind_inp
594	match_spec_model(len_index) = ind_mod
595	ENDIF
596
597	ENDIF
598
599	!
600	!-- Other Species
601	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) &
602	THEN
603	len_index = len_index + 1
604	match_spec_input(len_index) = ind_inp
605	match_spec_model(len_index) = ind_mod
606	ENDIF
607
608	END DO ! inp_ind
609
610	END DO ! inp_mod
611
612	!
613	!-- Error reporting (no matching)
614	ELSE
615
616	message_string = 'None of given Emission Species matches' // &
617	' model chemical species' // &
618	' Emission routine is not called'
619	CALL message( 'chem_emissions_matching', 'CM0440', 0, 0, 0, 6, 0 )
620
621	ENDIF
622
623	!
624	!-- Error reporting (no species)
625	ELSE
626
627	message_string = 'Array of Emission species not allocated: ' // &
628	' Either no emission species are provided as input or' // &
629	' no chemical species are used by PALM:' // &
630	' Emission routine is not called'
631	CALL message( 'chem_emissions_matching', 'CM0441', 0, 2, 0, 6, 0 )
632
633	ENDIF
634
635	!
636	!-- LOD 2 (PRE-PROCESSED mode)
637	CASE (2)
638
639	len_index = 0
640	len_index_voc = 0
641
642	IF ( nvar > 0 .AND. nspec_emis_inp > 0 ) THEN
643	!
644	!-- Cycle over model species
645	DO ind_mod = 1, nvar
646
647	!
648	!-- Cycle over input species
649	DO ind_inp = 1, nspec_emis_inp
650
651	!
652	!-- Check for VOC Species
653	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" ) THEN
654	DO ind_voc = 1, emt_att%nvoc
655	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) ) &
656	THEN
657	len_index = len_index + 1
658	len_index_voc = len_index_voc + 1
659	ENDIF
660	END DO
661	ENDIF
662
663	!
664	!-- Other Species
665	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) THEN
666	len_index = len_index + 1
667	ENDIF
668	ENDDO
669	ENDDO
670
671	!
672	!-- Allocate array for storing the indices of the matched species
673	IF ( len_index > 0 ) THEN
674
675	ALLOCATE ( match_spec_input(len_index) )
676
677	ALLOCATE ( match_spec_model(len_index) )
678
679	IF ( len_index_voc > 0 ) THEN
680	!
681	!-- Contains indices of the VOC model species
682	ALLOCATE( match_spec_voc_model(len_index_voc) )
683	!
684	!-- Contains the indices of different values of VOC composition of input variable
685	!-- VOC_composition
686	ALLOCATE( match_spec_voc_input(len_index_voc) )
687
688	ENDIF
689
690	!
691	!-- Pass the species indices to declared arrays
692	len_index = 0
693
694	!
695	!-- Cycle over model species
696	DO ind_mod = 1, nvar
697
698	!
699	!-- Cycle over Input species
700	DO ind_inp = 1, nspec_emis_inp
701
702	!
703	!-- VOCs
704	IF ( TRIM( emt_att%species_name(ind_inp) ) == "VOC" .AND. &
705	ALLOCATED( match_spec_voc_input ) ) THEN
706
707	DO ind_voc= 1, emt_att%nvoc
708	IF ( TRIM( emt_att%voc_name(ind_voc) ) == TRIM( spc_names(ind_mod) ) )&
709	THEN
710	len_index = len_index + 1
711	len_index_voc = len_index_voc + 1
712
713	match_spec_input(len_index) = ind_inp
714	match_spec_model(len_index) = ind_mod
715
716	match_spec_voc_input(len_index_voc) = ind_voc
717	match_spec_voc_model(len_index_voc) = ind_mod
718	ENDIF
719	END DO
720	ENDIF
721
722	!
723	!-- Other Species
724	IF ( TRIM( emt_att%species_name(ind_inp) ) == TRIM( spc_names(ind_mod) ) ) &
725	THEN
726	len_index = len_index + 1
727	match_spec_input(len_index) = ind_inp
728	match_spec_model(len_index) = ind_mod
729	ENDIF
730
731	END DO ! ind_inp
732	END DO ! ind_mod
733
734	ELSE ! if len_index_voc <= 0
735
736	!
737	!-- In case there are no species matching (just informational message)
738	message_string = 'Non of given emission species' // &
739	' matches' // &
740	' model chemical species:' // &
741	' Emission routine is not called'
742	CALL message( 'chem_emissions_matching', 'CM0438', 0, 0, 0, 6, 0 )
743	ENDIF
744
745	!
746	!-- Error check (no matching)
747	ELSE
748
749	!
750	!-- Either spc_names is zero or nspec_emis_inp is not allocated
751	message_string = 'Array of Emission species not allocated:' // &
752	' Either no emission species are provided as input or' // &
753	' no chemical species are used by PALM:' // &
754	' Emission routine is not called'
755	CALL message( 'chem_emissions_matching', 'CM0439', 0, 2, 0, 6, 0 )
756
757	ENDIF
758
759	!
760	!-- If emission module is switched on but mode_emis is not specified or it is given the wrong name
761
762	!
763	!-- Error check (no species)
764	CASE DEFAULT
765
766	message_string = 'Emission Module switched ON, but' // &
767	' either no emission mode specified or incorrectly given :' // &
768	' please, pass the correct value to the namelist parameter "mode_emis"'
769	CALL message( 'chem_emissions_matching', 'CM0445', 2, 2, 0, 6, 0 )
770
771	END SELECT
772
773	IF ( debug_output ) CALL debug_message( 'chem_emissions_match', 'end' )
774
775	END SUBROUTINE chem_emissions_match
776
777
778	!--------------------------------------------------------------------------------------------------!
779	! Description:
780	! ------------
781	!> Initialization:
782	!> Netcdf reading, arrays allocation and first calculation of cssws fluxes at timestep 0
783	!--------------------------------------------------------------------------------------------------!
784
785	SUBROUTINE chem_emissions_init
786
787	USE netcdf_data_input_mod, &
788	ONLY: chem_emis, chem_emis_att
789
790	IMPLICIT NONE
791
792	INTEGER(iwp) :: ispec !< running index
793
794	!
795	!-- Actions for initial runs
796	! IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
797	!-- ...
798	!
799	!
800	!-- Actions for restart runs
801	! ELSE
802	!-- ...
803	!
804	! ENDIF
805
806
807	IF ( debug_output ) CALL debug_message( 'chem_emissions_init', 'start' )
808
809	!
810	!-- Matching
811	CALL chem_emissions_match( chem_emis_att, n_matched_vars )
812
813	IF ( n_matched_vars == 0 ) THEN
814
815	emission_output_required = .FALSE.
816
817	ELSE
818
819	emission_output_required = .TRUE.
820
821
822	!
823	!-- Set molecule masses (in kg/mol)
824	ALLOCATE( chem_emis_att%xm(n_matched_vars) )
825
826	DO ispec = 1, n_matched_vars
827	SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
828	CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2
829	CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4
830	CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O
831	CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2
832	CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3
833	CASE ( 'CO' ); chem_emis_att%xm(ispec) = xm_C + xm_O
834	CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2
835	CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4
836	CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
837	CASE DEFAULT
838	chem_emis_att%xm(ispec) = 1.0_wp
839	END SELECT
840	ENDDO
841
842
843	!
844	!-- Get emissions for the first time step base on LOD (if defined) or emission mode
845	!-- (if no LOD defined)
846
847	!
848	!-- NOTE - I could use a combined if ( lod = xxx .or. mode = 'XXX' ) type of decision structure but
849	! I think it is much better to implement it this way (i.e., conditional on lod if it is
850	! defined, and mode if not) as we can easily take out the case structure for mode_emis
851	! later on.
852
853	IF ( emiss_lod < 0 ) THEN !-- no LOD defined (not likely)
854
855	SELECT CASE ( TRIM( mode_emis ) )
856
857	CASE ( 'PARAMETERIZED' ) ! LOD 0
858
859	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
860	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
861	ENDIF
862
863	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
864
865	CASE ( 'DEFAULT' ) ! LOD 1
866
867	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
868	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
869	ENDIF
870
871	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
872
873	CASE ( 'PRE-PROCESSED' ) ! LOD 2
874
875	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
876	!
877	!-- Note, at the moment emissions are considered only by surface fluxes rather than
878	!-- by volume sources. Therefore, no vertical dimension is required and is thus
879	!-- allocated with 1. Later when volume sources are considered, the vertical
880	!-- dimension will increase.
881	!ALLOCATE( emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,n_matched_vars) )
882	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
883	ENDIF
884
885	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
886
887	END SELECT
888
889	ELSE ! if LOD is defined
890
891	SELECT CASE ( emiss_lod )
892
893	CASE ( 0 ) ! parameterized mode
894
895	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
896	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
897	ENDIF
898
899	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars)
900
901	CASE ( 1 ) ! default mode
902
903	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
904	ALLOCATE( emis_distribution(1,nys:nyn,nxl:nxr,n_matched_vars) )
905	ENDIF
906
907	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
908
909	CASE ( 2 ) ! pre-processed mode
910
911	IF ( .NOT. ALLOCATED( emis_distribution) ) THEN
912	ALLOCATE( emis_distribution(nzb:nzt+1,nys:nyn,nxl:nxr,n_matched_vars) )
913	ENDIF
914
915	CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
916
917	END SELECT
918
919	ENDIF
920
921	!
922	! -- Initialize
923	emis_distribution = 0.0_wp
924
925	ENDIF
926
927	IF ( debug_output ) CALL debug_message( 'chem_emissions_init', 'end' )
928
929	END SUBROUTINE chem_emissions_init
930
931
932
933	!--------------------------------------------------------------------------------------------------!
934	! Description:
935	! ------------
936	!> Routine for Update of Emission values at each timestep.
937	!>
938	!> @todo Clarify the correct usage of index_dd, index_hh and index_mm. Consider renaming of these
939	!> variables.
940	!> @todo Clarify time used in emis_lod=2 mode. ATM, the used time seems strange.
941	!--------------------------------------------------------------------------------------------------!
942
943	SUBROUTINE chem_emissions_setup( emt_att, emt, n_matched_vars )
944
945	USE surface_mod, &
946	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
947
948	USE netcdf_data_input_mod, &
949	ONLY: street_type_f
950
951	USE arrays_3d, &
952	ONLY: hyp, pt
953
954	USE control_parameters, &
955	ONLY: time_since_reference_point
956
957	USE palm_date_time_mod, &
958	ONLY: days_per_week, get_date_time, hours_per_day, months_per_year, seconds_per_day
959
960	IMPLICIT NONE
961
962	INTEGER(iwp) :: day_of_month !< day of the month
963	INTEGER(iwp) :: day_of_week !< day of the week
964	INTEGER(iwp) :: day_of_year !< day of the year
965	INTEGER(iwp) :: days_since_reference_point !< days since reference point
966	INTEGER(iwp) :: i !< running index for grid in x-direction
967	INTEGER(iwp) :: i_pm_comp !< index for number of PM components
968	INTEGER(iwp) :: icat !< Index for number of categories
969	INTEGER(iwp) :: index_dd !< index day
970	INTEGER(iwp) :: index_hh !< index hour
971	INTEGER(iwp) :: index_mm !< index month
972	INTEGER(iwp) :: ispec !< index for number of species
973	INTEGER(iwp) :: ivoc !< Index for number of VOCs
974	INTEGER(iwp) :: hour_of_day !< hour of the day
975	INTEGER(iwp) :: j !< running index for grid in y-direction
976	INTEGER(iwp) :: k !< running index for grid in z-direction
977	INTEGER(iwp) :: m !< running index for horizontal surfaces
978	INTEGER(iwp) :: month_of_year !< month of the year
979
980	INTEGER,INTENT(IN) :: n_matched_vars !< Output of matching routine with number
981	!< of matched species
982
983	REAL(wp) :: time_utc_init !< second of day of initial date
984
985	TYPE(chem_emis_att_type), INTENT(INOUT) :: emt_att !< variable to store emission information
986
987	TYPE(chem_emis_val_type), INTENT(INOUT), ALLOCATABLE, DIMENSION(:) :: emt !< variable to store emission input values,
988	!< depending on the considered species
989	!
990	!-- CONVERSION FACTORS: TIME
991	REAL(wp), PARAMETER :: hour_per_year = 8760.0_wp !< number of hours in a year of 365 days
992	REAL(wp), PARAMETER :: s_per_hour = 3600.0_wp !< number of sec per hour (s)/(hour)
993	REAL(wp), PARAMETER :: s_per_day = 86400.0_wp !< number of sec per day (s)/(day)
994
995	REAL(wp), PARAMETER :: day_to_s = 1.0_wp/s_per_day !< conversion day -> sec
996	REAL(wp), PARAMETER :: hour_to_s = 1.0_wp/s_per_hour !< conversion hours -> sec
997	REAL(wp), PARAMETER :: year_to_s = 1.0_wp/(s_per_hour*hour_per_year) !< conversion year -> sec
998
999
1000	!
1001	!-- CONVERSION FACTORS: MASS
1002	REAL(wp), PARAMETER :: g_to_kg = 1.0E-03_wp !< Conversion from g to kg (kg/g)
1003	REAL(wp), PARAMETER :: miug_to_kg = 1.0E-09_wp !< Conversion from g to kg (kg/g)
1004	REAL(wp), PARAMETER :: tons_to_kg = 100.0_wp !< Conversion from tons to kg (kg/tons)
1005	!
1006	!-- CONVERSION FACTORS: PPM
1007	REAL(wp), PARAMETER :: ratio2ppm = 1.0E+06_wp
1008
1009	REAL(wp), DIMENSION(24) :: par_emis_time_factor !< time factors for the parameterized mode:
1010	!< fixed houlry profile for example day
1011	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: conv_to_ratio !< factor used for converting input
1012	!< to concentration ratio
1013	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: tmp_temp !< temporary variable for abs. temperature
1014
1015	REAL(wp), DIMENSION(:), ALLOCATABLE :: time_factor !< factor for time scaling of emissions
1016
1017	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: delta_emis !< incremental emission factor
1018	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: emis !< emission factor
1019
1020
1021	IF ( emission_output_required ) THEN
1022
1023	!
1024	!-- Set emis_dt to be used - since chemistry ODEs can be stiff, the option to solve them at every
1025	!-- RK substep is present to help improve stability should the need arises
1026
1027	IF ( call_chem_at_all_substeps ) THEN
1028
1029	dt_emis = dt_3d * weight_pres(intermediate_timestep_count)
1030
1031	ELSE
1032
1033	dt_emis = dt_3d
1034
1035	ENDIF
1036
1037	!
1038	!-- Conversion of units to the ones employed in PALM
1039	!-- In PARAMETERIZED mode no conversion is performed: in this case input units are fixed
1040	IF ( TRIM( mode_emis ) == "DEFAULT" .OR. TRIM( mode_emis ) == "PRE-PROCESSED" ) THEN
1041
1042	SELECT CASE ( TRIM( emt_att%units ) )
1043
1044	CASE ( 'kg/m2/s', 'KG/M2/S' ); conversion_factor = 1.0_wp ! kg
1045	CASE ( 'kg/m2/hour', 'KG/M2/HOUR' ); conversion_factor = hour_to_s
1046	CASE ( 'kg/m2/day', 'KG/M2/DAY' ); conversion_factor = day_to_s
1047	CASE ( 'kg/m2/year', 'KG/M2/YEAR' ); conversion_factor = year_to_s
1048
1049	CASE ( 'ton/m2/s', 'TON/M2/S' ); conversion_factor = tons_to_kg ! tonnes
1050	CASE ( 'ton/m2/hour', 'TON/M2/HOUR' ); conversion_factor = tons_to_kg*hour_to_s
1051	CASE ( 'ton/m2/year', 'TON/M2/YEAR' ); conversion_factor = tons_to_kg*year_to_s
1052
1053	CASE ( 'g/m2/s', 'G/M2/S' ); conversion_factor = g_to_kg ! grams
1054	CASE ( 'g/m2/hour', 'G/M2/HOUR' ); conversion_factor = g_to_kg*hour_to_s
1055	CASE ( 'g/m2/year', 'G/M2/YEAR' ); conversion_factor = g_to_kg*year_to_s
1056
1057	CASE ( 'micrograms/m2/s', 'MICROGRAMS/M2/S' ); conversion_factor = miug_to_kg ! ug
1058	CASE ( 'micrograms/m2/hour', 'MICROGRAMS/M2/HOUR' ); conversion_factor = miug_to_kg*hour_to_s
1059	CASE ( 'micrograms/m2/year', 'MICROGRAMS/M2/YEAR' ); conversion_factor = miug_to_kg*year_to_s
1060
1061	!
1062	!-- Error check (need units)
1063	CASE DEFAULT
1064	message_string = 'The Units of the provided emission input' // &
1065	' are not the ones required by PALM-4U: please check ' // &
1066	' emission module documentation.'
1067	CALL message( 'chem_emissions_setup', 'CM0446', 2, 2, 0, 6, 0 )
1068
1069	END SELECT
1070
1071	ENDIF
1072
1073	!
1074	!-- Conversion factor to convert kg/m**2/s to ppm/s
1075	DO i = nxl, nxr
1076	DO j = nys, nyn
1077
1078	!
1079	!-- Derive Temperature from Potential Temperature
1080	tmp_temp(nzb:nzt+1,j,i) = pt(nzb:nzt+1,j,i) * ( hyp(nzb:nzt+1) / p_0 )**rd_d_cp
1081
1082	!
1083	!-- We need to pass to cssws <- (ppm/s) * dz
1084	!-- Input is Nmole/(m^2*s)
1085	!-- To go to ppmdz multiply the input by (m2/N)dz
1086	!-- (m*2/N)dz == V/N
1087	!-- V/N = RT/P
1088	conv_to_ratio(nzb:nzt+1,j,i) = rgas_univ * & ! J K-1 mol-1
1089	tmp_temp(nzb:nzt+1,j,i) / & ! K
1090	hyp(nzb:nzt+1) ! Pa
1091
1092	! (ecc) for reference
1093	! m*3/Nmole (J/mol)K^-1 K Pa
1094	! conv_to_ratio(nzb:nzt+1,j,i) = ( (Rgas * tmp_temp(nzb:nzt+1,j,i)) / ((hyp(nzb:nzt+1))) )
1095
1096	ENDDO
1097	ENDDO
1098
1099
1100	! (ecc) moved initialization immediately after allocation
1101	!
1102	!-- Initialize
1103
1104	! emis_distribution(:,nys:nyn,nxl:nxr,:) = 0.0_wp
1105
1106
1107	!
1108	!-- LOD 2 (PRE-PROCESSED MODE)
1109
1110	IF ( emiss_lod == 2 ) THEN
1111
1112	! for reference (ecc)
1113	! IF ( TRIM( mode_emis ) == "PRE-PROCESSED" ) THEN
1114
1115	!
1116	!-- Update time indices
1117	CALL get_date_time( 0.0_wp, second_of_day=time_utc_init )
1118	CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ), hour=hour_of_day )
1119
1120	days_since_reference_point = INT( FLOOR( ( time_utc_init + &
1121	MAX( 0.0_wp, time_since_reference_point ) ) &
1122	/ seconds_per_day ) )
1123
1124	index_hh = days_since_reference_point * hours_per_day + hour_of_day
1125
1126	!
1127	!-- LOD 1 (DEFAULT MODE)
1128	ELSEIF ( emiss_lod == 1 ) THEN
1129
1130	! for reference (ecc)
1131	! ELSEIF ( TRIM( mode_emis ) == "DEFAULT" ) THEN
1132
1133	!
1134	!-- Allocate array where to store temporary emission values
1135	IF ( .NOT. ALLOCATED(emis) ) ALLOCATE( emis(nys:nyn,nxl:nxr) )
1136
1137	!
1138	!-- Allocate time factor per category
1139	ALLOCATE( time_factor(emt_att%ncat) )
1140
1141	!
1142	!-- Read-in hourly emission time factor
1143	IF ( TRIM( time_fac_type ) == "HOUR" ) THEN
1144
1145	!
1146	!-- Update time indices
1147	CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), &
1148	day_of_year=day_of_year, hour=hour_of_day )
1149	index_hh = ( day_of_year - 1_iwp ) * hour_of_day
1150
1151	!
1152	!-- Check if the index is less or equal to the temporal dimension of HOURLY emission files
1153	IF ( index_hh <= SIZE( emt_att%hourly_emis_time_factor(1,:) ) ) THEN
1154
1155	!
1156	!-- Read-in the correspondant time factor
1157	time_factor(:) = emt_att%hourly_emis_time_factor(:,index_hh+1)
1158
1159	!
1160	!-- Error check (time out of range)
1161	ELSE
1162
1163	message_string = 'The "HOUR" time-factors in the DEFAULT mode ' // &
1164	' are not provided for each hour of the total simulation time'
1165	CALL message( 'chem_emissions_setup', 'CM0448', 2, 2, 0, 6, 0 )
1166
1167	ENDIF
1168
1169	!
1170	!-- Read-in MDH emissions time factors
1171	ELSEIF ( TRIM( time_fac_type ) == "MDH" ) THEN
1172
1173	!
1174	!-- Update time indices
1175	CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), &
1176	month = month_of_year, &
1177	day = day_of_month, &
1178	hour = hour_of_day, &
1179	day_of_week = day_of_week &
1180	)
1181	index_mm = month_of_year
1182	index_dd = months_per_year + day_of_week
1183	SELECT CASE( TRIM( daytype_mdh ) )
1184
1185	CASE ("workday")
1186	index_hh = months_per_year + days_per_week + hour_of_day
1187
1188	CASE ("weekend")
1189	index_hh = months_per_year + days_per_week + hours_per_day + hour_of_day
1190
1191	CASE ("holiday")
1192	index_hh = months_per_year + days_per_week + 2*hours_per_day + hour_of_day
1193
1194	END SELECT
1195
1196	!
1197	!-- Check if the index is less or equal to the temporal dimension of MDH emission files
1198	IF ( ( index_hh + index_dd + index_mm) <= SIZE( emt_att%mdh_emis_time_factor(1,:) ) )&
1199	THEN
1200	!
1201	!-- Read in corresponding time factor
1202	time_factor(:) = emt_att%mdh_emis_time_factor(:,index_mm) * &
1203	emt_att%mdh_emis_time_factor(:,index_dd) * &
1204	emt_att%mdh_emis_time_factor(:,index_hh+1)
1205
1206	!
1207	!-- Error check (MDH time factor not provided)
1208	ELSE
1209
1210	message_string = 'The "MDH" time-factors in the DEFAULT mode ' // &
1211	' are not provided for each hour/day/month of the total simulation time'
1212	CALL message( 'chem_emissions_setup', 'CM0449', 2, 2, 0, 6, 0 )
1213
1214	ENDIF
1215
1216	!
1217	!-- Error check (no time factor defined)
1218	ELSE
1219
1220	message_string = 'In the DEFAULT mode the time factor' // &
1221	' has to be defined in the NAMELIST'
1222	CALL message( 'chem_emissions_setup', 'CM0450', 2, 2, 0, 6, 0 )
1223
1224	ENDIF
1225
1226	!
1227	!-- PARAMETERIZED MODE
1228	ELSEIF ( emiss_lod == 0 ) THEN
1229
1230
1231	! for reference (ecc)
1232	! ELSEIF ( TRIM( mode_emis ) == "PARAMETERIZED" ) THEN
1233
1234	!
1235	!-- Assign constant values of time factors, diurnal profile for traffic sector
1236	par_emis_time_factor(:) = (/ 0.009, 0.004, 0.004, 0.009, 0.029, 0.039, &
1237	0.056, 0.053, 0.051, 0.051, 0.052, 0.055, &
1238	0.059, 0.061, 0.064, 0.067, 0.069, 0.069, &
1239	0.049, 0.039, 0.039, 0.029, 0.024, 0.019 /)
1240
1241	IF ( .NOT. ALLOCATED( time_factor ) ) ALLOCATE( time_factor(1) )
1242
1243	!
1244	!-- Get time-factor for specific hour
1245	CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), hour = hour_of_day )
1246
1247	index_hh = hour_of_day
1248	time_factor(1) = par_emis_time_factor(index_hh+1)
1249
1250	ENDIF ! emiss_lod
1251
1252
1253	!
1254	!-- Emission distribution calculation
1255
1256	!
1257	!-- LOD 0 (PARAMETERIZED mode)
1258	IF ( emiss_lod == 0 ) THEN
1259
1260	! for reference (ecc)
1261	! IF ( TRIM( mode_emis ) == "PARAMETERIZED" ) THEN
1262
1263	DO ispec = 1, n_matched_vars
1264
1265	!
1266	!-- Units are micromoles/m*2day (or kilograms/m*2day for PMs)
1267	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = surface_csflux( match_spec_input(ispec) )&
1268	* time_factor(1) * hour_to_s
1269
1270	ENDDO
1271
1272
1273	!
1274	!-- LOD 1 (DEFAULT mode)
1275	ELSEIF ( emiss_lod == 1 ) THEN
1276
1277	! for referene (ecc)
1278	! ELSEIF ( TRIM( mode_emis ) == "DEFAULT" ) THEN
1279
1280	!
1281	!-- Allocate array for the emission value corresponding to a specific category and time factor
1282	ALLOCATE (delta_emis(nys:nyn,nxl:nxr))
1283
1284	!
1285	!-- Cycle over categories
1286	DO icat = 1, emt_att%ncat
1287
1288	!
1289	!-- Cycle over Species: n_matched_vars represents the number of species in common between
1290	!-- the emission input data and the chemistry mechanism used
1291	DO ispec = 1, n_matched_vars
1292
1293	emis(nys:nyn,nxl:nxr) = emt( match_spec_input(ispec) )% &
1294	default_emission_data(icat,nys+1:nyn+1,nxl+1:nxr+1)
1295
1296	!
1297	!-- NO
1298	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "NO" ) THEN
1299
1300	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
1301	time_factor(icat) * &
1302	emt_att%nox_comp(icat,1) * &
1303	conversion_factor * hours_per_day
1304
1305	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1306	emis_distribution(1,nys:nyn,nxl:nxr,ispec)&
1307	+ delta_emis(nys:nyn,nxl:nxr)
1308	!
1309	!-- NO2
1310	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "NO2" ) THEN
1311
1312	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
1313	time_factor(icat) * &
1314	emt_att%nox_comp(icat,2) * &
1315	conversion_factor * hours_per_day
1316
1317	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1318	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
1319	+ delta_emis(nys:nyn,nxl:nxr)
1320
1321	!
1322	!-- SO2
1323	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "SO2" ) THEN
1324
1325	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
1326	time_factor(icat) * &
1327	emt_att%sox_comp(icat,1) * &
1328	conversion_factor * hours_per_day
1329
1330	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1331	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
1332	+ delta_emis(nys:nyn,nxl:nxr)
1333
1334	!
1335	!-- SO4
1336	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "SO4" ) THEN
1337
1338
1339	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
1340	time_factor(icat) * &
1341	emt_att%sox_comp(icat,2) * &
1342	conversion_factor * hours_per_day
1343
1344	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1345	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
1346	+ delta_emis(nys:nyn,nxl:nxr)
1347
1348
1349	!
1350	!-- PM1
1351	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" ) THEN
1352
1353	DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,1) ) ! cycle through components
1354
1355	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
1356	time_factor(icat) * &
1357	emt_att%pm_comp(icat,i_pm_comp,1) * &
1358	conversion_factor * hours_per_day
1359
1360	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1361	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
1362	+ delta_emis(nys:nyn,nxl:nxr)
1363
1364	ENDDO
1365
1366	!
1367	!-- PM2.5
1368	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" ) THEN
1369
1370	DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,2) ) ! cycle through components
1371
1372	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
1373	time_factor(icat) * &
1374	emt_att%pm_comp(icat,i_pm_comp,2) * &
1375	conversion_factor * hours_per_day
1376
1377	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1378	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
1379	+ delta_emis(nys:nyn,nxl:nxr)
1380
1381	ENDDO
1382
1383	!
1384	!-- PM10
1385	ELSEIF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
1386
1387	DO i_pm_comp = 1, SIZE( emt_att%pm_comp(1,:,3) ) ! cycle through components
1388
1389	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * & ! kg/m2/s
1390	time_factor(icat) * &
1391	emt_att%pm_comp(icat,i_pm_comp,3) * &
1392	conversion_factor * hours_per_day
1393
1394	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1395	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
1396	+ delta_emis(nys:nyn,nxl:nxr)
1397
1398	ENDDO
1399
1400	!
1401	!-- VOCs
1402	ELSEIF ( SIZE( match_spec_voc_input ) > 0 ) THEN
1403
1404	DO ivoc = 1, SIZE( match_spec_voc_input ) ! cycle through components
1405
1406	IF ( TRIM( spc_names(match_spec_model(ispec) ) ) == &
1407	TRIM( emt_att%voc_name(ivoc) ) ) THEN
1408
1409	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * &
1410	time_factor(icat) * &
1411	emt_att%voc_comp(icat,match_spec_voc_input(ivoc)) * &
1412	conversion_factor * hours_per_day
1413
1414	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1415	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
1416	+ delta_emis(nys:nyn,nxl:nxr)
1417
1418	ENDIF
1419
1420	ENDDO
1421
1422	!
1423	!-- Any other species
1424	ELSE
1425
1426	delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr) * time_factor(icat) * &
1427	conversion_factor * hours_per_day
1428
1429	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1430	emis_distribution(1,nys:nyn,nxl:nxr,ispec) &
1431	+ delta_emis(nys:nyn,nxl:nxr)
1432
1433	ENDIF ! TRIM spc_names
1434
1435	emis = 0
1436
1437	ENDDO
1438
1439	delta_emis = 0
1440
1441	ENDDO
1442
1443	!
1444	!-- LOD 2 (PRE-PROCESSED mode)
1445	ELSEIF ( emiss_lod == 2 ) THEN
1446
1447	! for reference (ecc)
1448	! ELSEIF ( TRIM( mode_emis ) == "PRE-PROCESSED" ) THEN
1449
1450	!
1451	!-- Cycle over species: n_matched_vars represents the number of species in common between the
1452	!-- emission input data and the chemistry mechanism used
1453	DO ispec = 1, n_matched_vars
1454
1455	! (ecc)
1456	emis_distribution(1,nys:nyn,nxl:nxr,ispec) = emt(match_spec_input(ispec))% &
1457	preproc_emission_data(index_hh+1,1,nys+1:nyn+1,nxl+1:nxr+1) &
1458	* conversion_factor
1459
1460
1461	! emis_distribution(1,nys:nyn,nxl:nxr,ispec) = &
1462	! emt(match_spec_input(ispec))% &
1463	! preproc_emission_data(index_hh,1,:,:) * &
1464	! conversion_factor
1465	ENDDO
1466
1467	ENDIF ! emiss_lod
1468
1469
1470	!
1471	!-- Cycle to transform x,y coordinates to the one of surface_mod and to assign emission values to
1472	!-- cssws
1473
1474	!
1475	!-- LOD 0 (PARAMETERIZED mode)
1476	!-- Units of inputs are micromoles/m2/s
1477	IF ( emiss_lod == 0 ) THEN
1478	! for reference (ecc)
1479	! IF ( TRIM( mode_emis ) == "PARAMETERIZED" ) THEN
1480
1481	IF (street_type_f%from_file) THEN
1482
1483	!
1484	!-- Streets are lsm surfaces, hence, no usm surface treatment required.
1485	!-- However, urban surface may be initialized via default initialization in surface_mod, e.g. at
1486	!-- horizontal urban walls that are at k == 0 (building is lower than the first grid point). Hence,
1487	!-- in order to have only emissions at streets, set the surfaces emissions to zero at urban walls.
1488	IF ( surf_usm_h%ns >=1 ) surf_usm_h%cssws = 0.0_wp
1489
1490	!
1491	!-- Treat land-surfaces.
1492	DO m = 1, surf_lsm_h%ns
1493
1494	i = surf_lsm_h%i(m)
1495	j = surf_lsm_h%j(m)
1496	k = surf_lsm_h%k(m)
1497
1498	!
1499	!-- Set everything to zero then reassign according to street type
1500	surf_lsm_h%cssws(:,m) = 0.0_wp
1501
1502	IF ( street_type_f%var(j,i) >= main_street_id .AND. &
1503	street_type_f%var(j,i) < max_street_id ) THEN
1504
1505	!
1506	!-- Cycle over matched species
1507	DO ispec = 1, n_matched_vars
1508
1509	!
1510	!-- PMs are already in kilograms
1511	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
1512	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
1513	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
1514
1515	!
1516	!-- kg/(m^2s) kg/m^3
1517	surf_lsm_h%cssws(match_spec_model(ispec),m) = &
1518	emiss_factor_main(match_spec_input(ispec)) * &
1519	emis_distribution(1,j,i,ispec) * & ! kg/(m^2*s)
1520	rho_air(k) ! kg/m^3
1521
1522	!
1523	!-- Other Species
1524	!-- Inputs are micromoles
1525	ELSE
1526
1527	!
1528	!-- ppm/s m kg/m^3
1529	surf_lsm_h%cssws(match_spec_model(ispec),m) = &
1530	emiss_factor_main( match_spec_input(ispec) ) * &
1531	emis_distribution(1,j,i,ispec) * & ! micromoles/(m^2*s)
1532	conv_to_ratio(k,j,i) * & ! m^3/Nmole
1533	rho_air(k) ! kg/m^3
1534
1535	ENDIF
1536
1537	ENDDO ! ispec
1538
1539
1540	ELSEIF ( street_type_f%var(j,i) >= side_street_id .AND. &
1541	street_type_f%var(j,i) < main_street_id ) THEN
1542
1543	!
1544	!-- Cycle over matched species
1545	DO ispec = 1, n_matched_vars
1546
1547	!
1548	!-- PMs are already in kilograms
1549	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
1550	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
1551	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
1552
1553	!
1554	!-- kg/(m^2s) kg/m^3
1555	surf_lsm_h%cssws(match_spec_model(ispec),m) = &
1556	emiss_factor_side( match_spec_input(ispec) ) * &
1557	emis_distribution(1,j,i,ispec) * & ! kg/(m^2*s)
1558	rho_air(k) ! kg/m^3
1559	!
1560	!-- Other species
1561	!-- Inputs are micromoles
1562	ELSE
1563
1564	!
1565	!-- ppm/s m kg/m^3
1566	surf_lsm_h%cssws(match_spec_model(ispec),m) = &
1567	emiss_factor_side( match_spec_input(ispec) ) * &
1568	emis_distribution(1,j,i,ispec) * & ! micromoles/(m^2*s)
1569	conv_to_ratio(k,j,i) * & ! m^3/Nmole
1570	rho_air(k) ! kg/m^3
1571
1572	ENDIF
1573
1574	ENDDO ! ispec
1575
1576	!
1577	!-- If no street type is defined, then assign zero emission to all the species
1578	! (ecc) moved to front (for reference)
1579	! ELSE
1580	!
1581	! surf_lsm_h%cssws(:,m) = 0.0_wp
1582
1583	ENDIF ! street type
1584
1585	ENDDO ! m
1586
1587	ENDIF ! street_type_f%from_file
1588
1589
1590	!
1591	!-- LOD 1 (DEFAULT) and LOD 2 (PRE-PROCESSED)
1592	ELSE
1593
1594
1595	DO ispec = 1, n_matched_vars
1596
1597	!
1598	!-- Default surfaces
1599	DO m = 1, surf_def_h(0)%ns
1600
1601	i = surf_def_h(0)%i(m)
1602	j = surf_def_h(0)%j(m)
1603
1604	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
1605
1606	!
1607	!-- PMs
1608	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
1609	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
1610	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
1611
1612	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
1613	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
1614	rho_air(nzb) ! kg/m^3
1615
1616	ELSE
1617
1618	!
1619	!-- VOCs
1620	IF ( len_index_voc > 0 .AND. &
1621	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
1622
1623	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
1624	emis_distribution(1,j,i,ispec) * & ! mole/m2/s
1625	conv_to_ratio(nzb,j,i) * & ! m^3/mole
1626	ratio2ppm * & ! ppm
1627	rho_air(nzb) ! kg/m^3
1628
1629
1630	!
1631	!-- Other species
1632	ELSE
1633
1634	surf_def_h(0)%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
1635	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
1636	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
1637	conv_to_ratio(nzb,j,i) * & ! m^3/mole
1638	ratio2ppm * & ! ppm
1639	rho_air(nzb) ! kg/m^3
1640
1641	ENDIF ! VOC
1642
1643	ENDIF ! PM
1644
1645	ENDIF ! emis_distribution > 0
1646
1647	ENDDO ! m
1648
1649	!
1650	!-- LSM surfaces
1651	DO m = 1, surf_lsm_h%ns
1652
1653	i = surf_lsm_h%i(m)
1654	j = surf_lsm_h%j(m)
1655	k = surf_lsm_h%k(m)
1656
1657	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
1658
1659	!
1660	!-- PMs
1661	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
1662	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
1663	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
1664
1665	surf_lsm_h%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
1666	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
1667	rho_air(k) ! kg/m^3
1668
1669	ELSE
1670
1671	!
1672	!-- VOCs
1673	IF ( len_index_voc > 0 .AND. &
1674	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
1675
1676	surf_lsm_h%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
1677	emis_distribution(1,j,i,ispec) * & ! mole/m2/s
1678	conv_to_ratio(k,j,i) * & ! m^3/mole
1679	ratio2ppm * & ! ppm
1680	rho_air(k) ! kg/m^3
1681
1682	!
1683	!-- Other species
1684	ELSE
1685
1686	surf_lsm_h%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
1687	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
1688	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
1689	conv_to_ratio(k,j,i) * & ! m^3/mole
1690	ratio2ppm * & ! ppm
1691	rho_air(k) ! kg/m^3
1692
1693	ENDIF ! VOC
1694
1695	ENDIF ! PM
1696
1697	ENDIF ! emis_distribution
1698
1699	ENDDO ! m
1700
1701
1702
1703	!
1704	!-- USM surfaces
1705	DO m = 1, surf_usm_h%ns
1706
1707	i = surf_usm_h%i(m)
1708	j = surf_usm_h%j(m)
1709	k = surf_usm_h%k(m)
1710
1711	IF ( emis_distribution(1,j,i,ispec) > 0.0_wp ) THEN
1712
1713	!
1714	!-- PMs
1715	IF ( TRIM( spc_names( match_spec_model(ispec) ) ) == "PM1" .OR. &
1716	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM25" .OR. &
1717	TRIM( spc_names( match_spec_model(ispec) ) ) == "PM10" ) THEN
1718
1719	surf_usm_h%cssws(match_spec_model(ispec),m) = & ! kg/m2/s * kg/m3
1720	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
1721	rho_air(k) ! kg/m^3
1722
1723	ELSE
1724
1725	!
1726	!-- VOCs
1727	IF ( len_index_voc > 0 .AND. &
1728	emt_att%species_name(match_spec_input(ispec)) == "VOC" ) THEN
1729
1730	surf_usm_h%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
1731	emis_distribution(1,j,i,ispec) * & ! m2/s
1732	conv_to_ratio(k,j,i) * & ! m^3/mole
1733	ratio2ppm * & ! ppm
1734	rho_air(k) ! kg/m^3
1735
1736	!
1737	!-- Other species
1738	ELSE
1739
1740	surf_usm_h%cssws(match_spec_model(ispec),m) = & ! ppm/s * m * kg/m3
1741	emis_distribution(1,j,i,ispec) * & ! kg/m2/s
1742	( 1.0_wp / emt_att%xm(ispec) ) * & ! mole/kg
1743	conv_to_ratio(k,j,i) * & ! m^3/mole
1744	ratio2ppm * & ! ppm
1745	rho_air(k) ! kg/m^3
1746
1747
1748	ENDIF ! VOC
1749
1750	ENDIF ! PM
1751
1752	ENDIF ! emis_distribution
1753
1754	ENDDO ! m
1755
1756	ENDDO
1757
1758	ENDIF
1759
1760	!
1761	!-- Deallocate time_factor in case of DEFAULT mode)
1762	IF ( ALLOCATED( time_factor ) ) DEALLOCATE( time_factor )
1763
1764	ENDIF
1765
1766	END SUBROUTINE chem_emissions_setup
1767
1768
1769	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1770	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1771	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1772	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1773	!!
1774	!! 20200203 (ECC) - ON DEMAND EMISSION UPDATE MODE
1775	!!
1776	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1777	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1778	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1779	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1780
1781
1782	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1783	!!
1784	!! WRAPPER / INTERFACE FUNCTIONS
1785	!!
1786	!! NOTE - I find using an explicity wrapper provides much better flow control
1787	!! over an interface block
1788	!!
1789	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1790
1791	!
1792	!-- 20200203 (ECC)
1793	!
1794	!--------------------------------------------------------------------------------------------------!
1795	! Description:
1796	! ------------
1797	!> interface for initiation of emission arrays based on emission LOD
1798	!
1799	!--------------------------------------------------------------------------------------------------!
1800
1801	SUBROUTINE chem_emissions_header_init
1802
1803	IMPLICIT NONE
1804
1805	SELECT CASE ( emiss_lod )
1806	CASE ( 0 )
1807	! do nothing at the moment
1808	CASE ( 1 )
1809	! do nothing at the moment
1810	CASE ( 2 )
1811	CALL chem_emissions_header_init_lod2
1812	END SELECT
1813
1814	END SUBROUTINE chem_emissions_header_init
1815
1816
1817	!
1818	!-- 20200203 (ECC)
1819	!
1820	!--------------------------------------------------------------------------------------------------!
1821	! Description:
1822	! ------------
1823	!> interface for initiation of emission arrays based on emission LOD
1824	!--------------------------------------------------------------------------------------------------!
1825
1826	SUBROUTINE chem_emissions_update_on_demand
1827
1828	IMPLICIT NONE
1829
1830	SELECT CASE ( emiss_lod )
1831	CASE ( 0 )
1832	! do nothing at the moment
1833	CASE ( 1 )
1834	! do nothing at the moment
1835	CASE ( 2 )
1836	CALL chem_emissions_update_on_demand_lod2
1837	END SELECT
1838
1839	END SUBROUTINE ! chem_emisisons_update_on_demand
1840
1841
1842	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1843	!!
1844	!! SUBROUTINES SPECIFIC FOR LOD 2
1845	!!
1846	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1847
1848	!
1849	!-- 20200203 (ECC)
1850	!
1851	!--------------------------------------------------------------------------------------------------!
1852	! Description:
1853	! ------------
1854	!> Initiates header for emissions data attributes for LOD 2
1855	!--------------------------------------------------------------------------------------------------!
1856
1857	SUBROUTINE chem_emissions_header_init_lod2
1858
1859	USE control_parameters, &
1860	ONLY: coupling_char, message_string
1861
1862	USE netcdf_data_input_mod, &
1863	ONLY: chem_emis_att, close_input_file, get_attribute, get_dimension_length, get_variable, &
1864	open_read_file
1865
1866
1867	IMPLICIT NONE
1868
1869
1870	INTEGER(iwp) :: att_lod !< lod attribute in chemistry file
1871	INTEGER(iwp) :: ncid !< chemistry file netCDF handle
1872
1873	IF ( debug_output ) CALL debug_message( 'chem_emissions_header_init_lod2', 'start' )
1874
1875	!
1876	!-- Opens _chemistry input file and obtain header information
1877	CALL open_read_file ( TRIM( input_file_chem ) // TRIM( coupling_char ), ncid )
1878	!
1879	!-- Check if LOD in chemistry file matches LOD in namelist
1880	CALL get_attribute ( ncid, 'lod', att_lod, .TRUE. )
1881
1882	IF ( att_lod /= emiss_lod ) THEN
1883	message_string = '' ! to get around unused variable warning / error
1884	WRITE ( message_string, * ) 'LOD mismatch between namelist (emiss_lod) and', &
1885	CHAR( 10 ), ' ', 'chemistry input file (global attributes>lod)'
1886	CALL message( 'chem_emissions_header_init_lod2', 'CM0468', 1, 2, 0, 6, 0 )
1887	ENDIF
1888	!
1889	!-- Obtain unit conversion factor
1890	CALL get_attribute ( ncid, 'units', chem_emis_att%units, .FALSE., "emission_values" )
1891	conversion_factor = chem_emissions_convert_base_units ( chem_emis_att%units )
1892	!
1893	!-- Obtain header attributes
1894	CALL chem_emissions_init_species ( ncid )
1895	CALL chem_emissions_init_timestamps ( ncid )
1896	!
1897	!-- Done reading file
1898	CALL close_input_file (ncid)
1899
1900	!
1901	!-- Set previous timestamp index to something different to trigger a read event later on
1902	previous_timestamp_index = -1
1903
1904	IF ( debug_output ) CALL debug_message( 'chem_emissions_header_init_lod2', 'end' )
1905
1906	END SUBROUTINE chem_emissions_header_init_lod2
1907
1908	!
1909	!-- 20200203 (ECC)
1910	!
1911	!--------------------------------------------------------------------------------------------------!
1912	! Description:
1913	! ------------
1914	!> Reads emission data on demand for LOD2
1915	!--------------------------------------------------------------------------------------------------!
1916
1917	SUBROUTINE chem_emissions_update_on_demand_lod2
1918
1919	USE control_parameters, &
1920	ONLY: coupling_char, time_since_reference_point
1921
1922	USE netcdf_data_input_mod, &
1923	ONLY: chem_emis_att, close_input_file, get_variable, open_read_file
1924
1925	USE arrays_3d, &
1926	ONLY: hyp, pt
1927
1928	USE surface_mod, &
1929	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
1930
1931
1932	IMPLICIT NONE
1933
1934	CHARACTER(LEN=80) :: this_timestamp !< writes out timestamp
1935
1936	INTEGER(iwp) :: i,j,k,m !< generic counters
1937	INTEGER(iwp) :: kmatch !< index of matched species
1938	INTEGER(iwp) :: ncid !< netCDF file handle (chemistry file)
1939	INTEGER(iwp) :: time_index_location !< location of most recent timestamp
1940
1941	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: cssws_def_h !< dummy default surface array
1942	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: cssws_lsm_h !< dummy LSM surface array
1943	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: cssws_usm_h !< dummy USM surface array
1944	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: mass2mole !< conversion factor mass 2 molar (ppm) flux
1945
1946	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: emis_distrib !< surface emissions
1947
1948	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:) :: emissions_raw !< raw emissions data
1949
1950	IF ( debug_output ) CALL debug_message ( 'chem_emissions_update_on_demand_lod2', 'start' )
1951	!
1952	!-- Obtain current timestamp and locate index for most recent timestamp element
1953	!-- end subroutine (RETURN) if it is still the same index as the existing time index
1954
1955	this_timestamp = '' ! string must be initiated before using
1956	CALL get_date_time( time_since_reference_point, date_time_str=this_timestamp )
1957
1958	time_index_location = chem_emissions_locate_timestep &
1959	( this_timestamp, timestamps, 1, chem_emis_att%dt_emission )
1960
1961	IF ( time_index_location == previous_timestamp_index ) RETURN
1962
1963	!
1964	!-- Begin extract emission data for matched species from netCDF file
1965	previous_timestamp_index = time_index_location
1966
1967	ALLOCATE ( emis_distrib(n_matched_vars,nys:nyn,nxl:nxr) )
1968	emis_distrib = 0.0_wp
1969
1970	!
1971	!-- Open netCDF file and allocate temp memory
1972	CALL open_read_file( TRIM( input_file_chem ) // TRIM( coupling_char ), ncid )
1973	ALLOCATE( emissions_raw(1,1,nys:nyn,nxl:nxr,1) )
1974
1975	DO k = 1, n_matched_vars
1976	!
1977	!-- Get index for matching species
1978	kmatch = chem_emissions_locate_species( spc_names(match_spec_model(k)), &
1979	chem_emis_att%species_name )
1980	!
1981	!-- Extract variable as-is
1982	!-- Note C index notations for nx and ny due to MPI and reversed index dimension order for netCDF
1983	!-- Fortran API)
1984	emissions_raw = 0.0_wp
1985
1986	CALL get_variable ( ncid, 'emission_values', emissions_raw, &
1987	kmatch, nxl+1, nys+1, 1, time_index_location, &
1988	1, nxr-nxl+1, nyn-nys+1, 1, 1, .FALSE. )
1989	!
1990	!-- Transfer emission data
1991	DO j = nys,nyn
1992	DO i = nxl,nxr
1993	emis_distrib(k,j,i) = emissions_raw(1,1,j,i,1) * conversion_factor
1994	ENDDO
1995	ENDDO
1996
1997	ENDDO ! k = n_matched_vars
1998	!
1999	!-- netCDF handle and temp memory no longer needed
2000	DEALLOCATE( emissions_raw )
2001	CALL close_input_file( ncid )
2002	!
2003	!-- Set emis_dt since chemistry ODEs can be stiff, the option to solve them at every RK substep is
2004	!-- present to help improve stability should the need arise
2005	dt_emis = dt_3d
2006
2007	IF ( call_chem_at_all_substeps ) dt_emis = dt_emis * weight_pres(intermediate_timestep_count)
2008	!
2009	!-- Calculate conversion factor from mass flux to molar flux (mixing ratio)
2010	ALLOCATE ( mass2mole(nys:nyn,nxl:nxr) )
2011	mass2mole = 0.0_wp
2012
2013	DO i = nxl, nxr
2014	DO j = nys, nyn
2015	mass2mole(j,i) = mass_2_molar_flux ( hyp(nzb), pt(nzb,j,i) )
2016	ENDDO
2017	ENDDO
2018
2019	!
2020	!-- Calculate surface fluxes
2021	!-- NOTE - For some reason I can not pass surf_xxx%cssws as output argument into subroutine
2022	!-- assign_surface_flux ( ). The contents got mixed up once the subroutine is finished. I
2023	!-- don't know why and I don't have time to investigate. As workaround I declared dummy
2024	!-- variables and reassign them one by one (i.e., in a loop)
2025	!-- ECC 20200206
2026
2027	!
2028	!-- Allocate and initialize dummy surface fluxes
2029	ALLOCATE( cssws_def_h(n_matched_vars,surf_def_h(0)%ns) )
2030	cssws_def_h = 0.0_wp
2031
2032	ALLOCATE( cssws_lsm_h(n_matched_vars,surf_lsm_h%ns) )
2033	cssws_lsm_h = 0.0_wp
2034
2035	ALLOCATE( cssws_usm_h(n_matched_vars,surf_usm_h%ns) )
2036	cssws_usm_h = 0.0_wp
2037
2038	!
2039	!-- Assign and transfer emissions as surface fluxes
2040	CALL assign_surface_flux ( cssws_def_h, surf_def_h(0)%ns, &
2041	surf_def_h(0)%j, surf_def_h(0)%i, &
2042	emis_distrib, mass2mole )
2043
2044
2045	CALL assign_surface_flux ( cssws_lsm_h, surf_lsm_h%ns, &
2046	surf_lsm_h%j, surf_lsm_h%i, &
2047	emis_distrib, mass2mole )
2048
2049
2050	CALL assign_surface_flux ( cssws_usm_h, surf_usm_h%ns, &
2051	surf_usm_h%j, surf_usm_h%i, &
2052	emis_distrib, mass2mole )
2053
2054	DO k = 1, n_matched_vars
2055
2056	DO m = 1, surf_def_h(0)%ns
2057	surf_def_h(0)%cssws(k,m) = cssws_def_h(k,m)
2058	ENDDO
2059
2060	DO m = 1, surf_lsm_h%ns
2061	surf_lsm_h%cssws(k,m) = cssws_lsm_h(k,m)
2062	ENDDO
2063
2064	DO m = 1, surf_usm_h%ns
2065	surf_usm_h%cssws(k,m) = cssws_usm_h(k,m)
2066	ENDDO
2067
2068	ENDDO
2069
2070	!
2071	!-- Cleaning up
2072	DEALLOCATE( cssws_def_h )
2073	DEALLOCATE( cssws_lsm_h )
2074	DEALLOCATE( cssws_usm_h )
2075
2076	DEALLOCATE ( emis_distrib )
2077	DEALLOCATE ( mass2mole )
2078
2079	IF ( debug_output ) CALL debug_message ( 'chem_emissions_update_on_demand_lod2', 'end' )
2080
2081	END SUBROUTINE ! chem_emissions_update_on_demand_lod2
2082
2083
2084	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2085	!!
2086	!! AUXILIARY SUBROUTINES
2087	!!
2088	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2089
2090	!
2091	!-- 20200203 (ECC)
2092	!
2093	!--------------------------------------------------------------------------------------------------!
2094	! Description:
2095	! ------------
2096	!> Look for matched species between emissions attributes and selected chemical mechanisms and
2097	!> determine corresponding molecular weights
2098	!--------------------------------------------------------------------------------------------------!
2099
2100	SUBROUTINE chem_emissions_init_species ( ncid )
2101
2102	USE netcdf_data_input_mod, &
2103	ONLY: chem_emis_att, close_input_file, get_dimension_length, get_variable, open_read_file
2104
2105	IMPLICIT NONE
2106
2107	INTEGER(iwp) :: ispec !< generic counter 4 species
2108
2109	INTEGER(iwp), INTENT(IN) :: ncid !< netcdf file ID
2110
2111	IF ( debug_output ) CALL debug_message( 'chem_emissions_init_species', 'start' )
2112	!
2113	!- Assign species
2114	CALL get_dimension_length ( ncid, chem_emis_att%n_emiss_species, 'nspecies' )
2115	CALL get_variable ( ncid, 'emission_name', chem_emis_att%species_name, &
2116	chem_emis_att%n_emiss_species )
2117	!
2118	!- Backward compatibility for salsa_mod (ECC)
2119	chem_emis_att%nspec = chem_emis_att%n_emiss_species
2120	!
2121	!-- Get a list of matched species between emission_attributes and selected chemical mechanism
2122	emission_output_required = .FALSE.
2123	CALL chem_emissions_match( chem_emis_att, n_matched_vars )
2124
2125	!
2126	!-- If matched species found (at least 1),
2127	!-- allocate memory for emission attributes,
2128	!-- assign molecular masses [kg/mol],
2129	!-- see chemistry_model_mod.f90 for reference.
2130	IF ( n_matched_vars > 0 ) THEN
2131
2132	emission_output_required = .TRUE.
2133
2134	ALLOCATE( chem_emis_att%xm(n_matched_vars) )
2135
2136	DO ispec = 1, n_matched_vars
2137	chem_emis_att%xm(ispec) = 1.0_wp
2138	SELECT CASE ( TRIM( spc_names(match_spec_model(ispec)) ) )
2139	CASE ( 'SO2' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 2
2140	CASE ( 'SO4' ); chem_emis_att%xm(ispec) = xm_S + xm_O * 4
2141	CASE ( 'NO' ); chem_emis_att%xm(ispec) = xm_N + xm_O
2142	CASE ( 'NO2' ); chem_emis_att%xm(ispec) = xm_N + xm_O * 2
2143	CASE ( 'NH3' ); chem_emis_att%xm(ispec) = xm_N + xm_H * 3
2144	CASE ( 'CO' ); chem_emis_att%xm(ispec) = xm_C + xm_O
2145	CASE ( 'CO2' ); chem_emis_att%xm(ispec) = xm_C + xm_O * 2
2146	CASE ( 'CH4' ); chem_emis_att%xm(ispec) = xm_C + xm_H * 4
2147	CASE ( 'HNO3' ); chem_emis_att%xm(ispec) = xm_H + xm_N + xm_O*3
2148	END SELECT
2149	ENDDO
2150
2151	ENDIF ! IF ( n_matched_vars > 0 )
2152
2153	IF ( debug_output ) CALL debug_message( 'chem_emissions_init_species', 'end' )
2154
2155	END SUBROUTINE chem_emissions_init_species
2156
2157
2158	!
2159	!-- 20200203 (ECC)
2160	!
2161	!--------------------------------------------------------------------------------------------------!
2162	! Description:
2163	! ------------
2164	!> Extract timestamps from netCDF input
2165	!--------------------------------------------------------------------------------------------------!
2166
2167	SUBROUTINE chem_emissions_init_timestamps ( ncid )
2168
2169	USE control_parameters, &
2170	ONLY: message_string
2171
2172	USE netcdf_data_input_mod, &
2173	ONLY: chem_emis_att, close_input_file, get_dimension_length, get_variable, open_read_file
2174
2175	IMPLICIT NONE
2176
2177	INTEGER(iwp) :: fld_len !< string field length
2178	INTEGER(iwp) :: itime !< generic counter (4 species)
2179
2180	INTEGER(iwp), INTENT(IN) :: ncid !< netcdf file handle
2181
2182	IF ( debug_output ) CALL debug_message( 'chem_emissions_read_timestamps', 'start' )
2183	!
2184	!-- Import timestamps from netCDF input
2185	CALL get_dimension_length ( ncid, chem_emis_att%dt_emission, 'time' )
2186	CALL get_dimension_length ( ncid, fld_len, 'field_length' )
2187	CALL get_variable ( ncid, 'timestamp', timestamps, chem_emis_att%dt_emission, fld_len )
2188	!
2189	!-- Throw error at first instance of timestamps not in listed in chronological order.
2190	DO itime = 2,chem_emis_att%dt_emission
2191
2192	IF ( timestamps(itime-1) > timestamps(itime) ) THEN
2193
2194	WRITE( message_string, * ) &
2195	'input timestamps not in chronological order for', &
2196	CHAR( 10 ), ' ', &
2197	'index ', (itime-1), ' : ', TRIM( timestamps(itime-1) ), ' and', &
2198	CHAR( 10 ), ' ', &
2199	'index ', (itime), ' : ', TRIM( timestamps(itime) )
2200
2201	CALL message( 'chem_emissions_read_timestamps', 'CM0469', 1, 2, 0, 6, 0 )
2202
2203	ENDIF
2204
2205	ENDDO
2206
2207	IF ( debug_output ) CALL debug_message( 'chem_emissions_read_timestamps', 'end' )
2208
2209	END SUBROUTINE chem_emissions_init_timestamps
2210
2211
2212	!
2213	!-- 20200203 (ECC)
2214	!
2215	!--------------------------------------------------------------------------------------------------!
2216	! Description:
2217	! ------------
2218	!> Assign emissions as surface fluxes
2219	!
2220	!> NOTE: For arguments, I originally wanted to use unspecified dimensions, but I could not get
2221	!> this to work properly, hence the dimensioned array arguments.
2222	!--------------------------------------------------------------------------------------------------!
2223
2224	SUBROUTINE assign_surface_flux ( surf_array, nsurfs, surf_j, surf_i, emis_dist, conv_mole )
2225
2226	USE arrays_3d, &
2227	ONLY: rho_air
2228
2229	USE netcdf_data_input_mod, &
2230	ONLY: chem_emis_att
2231
2232	USE surface_mod !< for surf_type
2233
2234	IMPLICIT NONE
2235	!
2236	!-- Input arguments
2237	INTEGER(iwp), INTENT(IN) :: nsurfs !< # surfaces in surf_array
2238
2239	INTEGER(iwp), DIMENSION(nsurfs), INTENT(IN) :: surf_i !< i indices 4 surf. elements
2240	INTEGER(iwp), DIMENSION(nsurfs), INTENT(IN) :: surf_j !< j indices 4 surf. elements
2241
2242	REAL(wp), DIMENSION(nys:nyn,nxl:nxr), INTENT(IN) :: conv_mole !< conv. 2 molar flux
2243	REAL(wp), DIMENSION(n_matched_vars,nys:nyn,nxl:nxr), INTENT(IN) :: emis_dist !< surf. emissions
2244
2245	REAL(wp), DIMENSION(n_matched_vars,nsurfs), INTENT(INOUT) :: surf_array !< surface listing
2246
2247	!
2248	!-- Parameters (magic numbers)
2249	CHARACTER(LEN=2), PARAMETER :: sp_PM = 'PM' !< id string for all PMs
2250	CHARACTER(LEN=3), PARAMETER :: sp_VOC = 'VOC' !< id string for VOC
2251
2252	REAL(wp), PARAMETER :: mol2ppm = 1.0E+06_wp !< conversion from mole 2 ppm
2253	!
2254	!-- Local variables
2255	CHARACTER(LEN=80) :: this_species_name !< matched species name
2256
2257	INTEGER(iwp) :: i,j,k,m !< generic counters
2258
2259	REAL(wp) :: flux_conv_factor !< conversion factor
2260
2261	IF ( debug_output ) CALL debug_message( 'chem_emissions_header_init_lod2', 'start' )
2262
2263	DO k = 1, n_matched_vars
2264
2265	this_species_name = spc_names(k) !< species already matched
2266
2267	DO m = 1, nsurfs
2268
2269	j = surf_j(m) ! get surface coordinates
2270	i = surf_i(m)
2271
2272	!
2273	!-- Calculate conversion factor depending on emission species type
2274	flux_conv_factor = rho_air(nzb)
2275	!
2276	!-- Account for conversion to different types of emisison species
2277	IF ( TRIM( this_species_name( 1:LEN( sp_PM ) ) ) == sp_PM ) THEN
2278
2279	! do nothing (use mass flux directly)
2280
2281	ELSE IF ( TRIM( this_species_name( 1:LEN( sp_VOC ) ) ) == sp_VOC ) THEN
2282
2283	flux_conv_factor = flux_conv_factor * conv_mole(j,i) * mol2ppm
2284
2285	ELSE
2286
2287	flux_conv_factor = flux_conv_factor * conv_mole(j,i) * mol2ppm / chem_emis_att%xm(k)
2288
2289	ENDIF
2290	!
2291	!-- Finally assign surface flux
2292	surf_array(k,m) = emis_dist(k,j,i) * flux_conv_factor
2293
2294	ENDDO ! m = 1, nsurfs
2295
2296	ENDDO ! k = 1, n_matched_vars
2297
2298
2299	IF ( debug_output ) CALL debug_message( 'chem_emissions_header_init_lod2', 'end' )
2300
2301	END SUBROUTINE assign_surface_flux
2302
2303
2304	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2305	!!
2306	!! AUXILIARY FUNCTIONS
2307	!!
2308	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2309
2310	!
2311	!-- 20200203 (ECC)
2312	!
2313	!--------------------------------------------------------------------------------------------------!
2314	! Description:
2315	! ------------
2316	!> Given incoming flux units ( mass / area / time ) provide single-valued onversion factor to
2317	!> ( kg / m2 / s )
2318	!--------------------------------------------------------------------------------------------------!
2319
2320	FUNCTION chem_emissions_convert_base_units ( units_in ) RESULT ( conv_factor )
2321
2322	IMPLICIT NONE
2323	!
2324	!-- Function arguments
2325	CHARACTER(LEN=*), INTENT(IN) :: units_in !< incoming units (ie emt_att%units)
2326
2327	REAL(wp) :: conv_factor !< convertion factor
2328
2329	!
2330	!-- Parameters (magic numbers)
2331	INTEGER(iwp), PARAMETER :: up2lo = 32 !< convert letter to lower case
2332	!
2333	!-- Base unit conversion factors (should be self-explanatory)
2334	REAL(wp), PARAMETER :: hour_per_year = 8760.0_wp
2335	REAL(wp), PARAMETER :: g_to_kg = 1.0E-03_wp
2336	REAL(wp), PARAMETER :: miug_to_kg = 1.0E-09_wp
2337	REAL(wp), PARAMETER :: s_per_hour = 3600.0_wp
2338	REAL(wp), PARAMETER :: s_per_day = 86400.0_wp
2339	REAL(wp), PARAMETER :: tons_to_kg = 100.0_wp
2340
2341	REAL(wp), PARAMETER :: day_to_s = 1.0_wp / s_per_day
2342	REAL(wp), PARAMETER :: hour_to_s = 1.0_wp / s_per_hour
2343	REAL(wp), PARAMETER :: year_to_s = 1.0_wp / s_per_hour / hour_per_year
2344
2345
2346	!
2347	!-- Local variables
2348	CHARACTER(LEN=LEN(units_in)) :: units_in_lo !< units in lower case
2349
2350	INTEGER(iwp) :: j,k !< generic counters
2351	INTEGER(iwp) :: str_len !< length of unit string
2352	!
2353	!-- Turn units string to lower case
2354	units_in_lo = ''
2355	str_len = LEN( TRIM( units_in ) )
2356
2357	DO k = 1,str_len
2358	j = IACHAR( units_in(k:k) )
2359	units_in_lo(k:k) = ACHAR( j )
2360	IF ( ( j >= IACHAR( "A" ) ) .AND. ( j <= IACHAR( "Z" ) ) ) &
2361	units_in_lo(k:k) = ACHAR( j + up2lo )
2362	ENDDO
2363
2364	conv_factor = 1.0_wp !< default value (kg/m2/s)
2365
2366	SELECT CASE ( TRIM( units_in_lo ) )
2367	CASE ( 'kg/m2/s' ); conv_factor = 1.0_wp
2368	CASE ( 'kg/m2/hour' ); conv_factor = hour_to_s
2369	CASE ( 'kg/m2/day' ); conv_factor = day_to_s
2370	CASE ( 'kg/m2/year' ); conv_factor = year_to_s
2371	CASE ( 'ton/m2/s' ); conv_factor = tons_to_kg
2372	CASE ( 'ton/m2/hour' ); conv_factor = tons_to_kg * hour_to_s
2373	CASE ( 'ton/m2/day' ); conv_factor = tons_to_kg * day_to_s
2374	CASE ( 'ton/m2/year' ); conv_factor = tons_to_kg * year_to_s
2375	CASE ( 'g/m2/s' ); conv_factor = g_to_kg
2376	CASE ( 'g/m2/hour' ); conv_factor = g_to_kg * hour_to_s
2377	CASE ( 'g/m2/day' ); conv_factor = g_to_kg * day_to_s
2378	CASE ( 'g/m2/year' ); conv_factor = g_to_kg * year_to_s
2379	CASE ( 'micrograms/m2/s' ); conv_factor = miug_to_kg
2380	CASE ( 'micrograms/m2/hour' ); conv_factor = miug_to_kg * hour_to_s
2381	CASE ( 'micrograms/m2/day' ); conv_factor = miug_to_kg * day_to_s
2382	CASE ( 'micrograms/m2/year' ); conv_factor = miug_to_kg * year_to_s
2383	CASE DEFAULT
2384	message_string = '' ! to get around unused variable warning / error
2385	WRITE ( message_string, * ) 'Specified emission units (', TRIM( units_in ), &
2386	') not recognized in PALM-4U'
2387	CALL message ( 'chem_emission_convert_units', 'CM0446', 2, 2, 0, 6, 0 )
2388	END SELECT
2389
2390	END FUNCTION chem_emissions_convert_base_units
2391
2392
2393	!
2394	!-- 20200203 (ECC)
2395	!
2396	!--------------------------------------------------------------------------------------------------!
2397	! Description:
2398	! ------------
2399	!> Calculates conversion factor from mass flux to ppm (molar flux)
2400	!--------------------------------------------------------------------------------------------------!
2401
2402	FUNCTION mass_2_molar_flux ( rhogh, theta ) RESULT ( conv_factor )
2403
2404	USE basic_constants_and_equations_mod, &
2405	ONLY: p_0, rd_d_cp, rgas_univ
2406
2407	IMPLICIT NONE
2408	!
2409	!-- Function arguments
2410	REAL(wp) :: conv_factor !< conversion factor
2411	REAL(wp), INTENT(IN) :: rhogh !< hydrostatic pressure
2412	REAL(wp), INTENT(IN) :: theta !< potential temperature
2413
2414	conv_factor = ( rgas_univ / rhogh ) * theta * ( ( rhogh / p_0 ) ** rd_d_cp )
2415
2416	END FUNCTION mass_2_molar_flux
2417
2418
2419	!
2420	!-- 20200203 (ECC)
2421	!
2422	!--------------------------------------------------------------------------------------------------!
2423	! Description:
2424	! ------------
2425	!> Given target sepecies locate index in species array
2426	!> returns 0 if none is found
2427	!--------------------------------------------------------------------------------------------------!
2428
2429	FUNCTION chem_emissions_locate_species ( this_species, species_array ) RESULT ( species_index )
2430
2431	IMPLICIT NONE
2432	!
2433	!-- Function arguments
2434	INTEGER(iwp) :: species_index !> index matching species
2435
2436	CHARACTER(LEN=25), INTENT(IN) :: species_array(:) !> array of species
2437	CHARACTER(LEN=*), INTENT(IN) :: this_species !> target species
2438	!
2439	!-- Local variables
2440	INTEGER(iwp) :: k !> generic counter
2441	INTEGER(iwp) :: n_species !> number of species in species_array
2442
2443	n_species = SIZE( species_array, 1 )
2444
2445	DO k = 1, n_species
2446	IF ( TRIM( species_array(k) ) == TRIM( this_species ) ) EXIT
2447	ENDDO
2448
2449	species_index = 0 !> assume no matching index is found
2450
2451	IF ( TRIM( species_array(k) ) == TRIM( this_species ) ) specieS_index = k
2452
2453	END FUNCTION chem_emissions_locate_species
2454
2455
2456	!
2457	!-- 20200203 (ECC)
2458	!
2459	!--------------------------------------------------------------------------------------------------!
2460	! Description:
2461	! ------------
2462	!> given target timestamp locate most recent timestep in timestamp array
2463	!> using bisection search (since array is sorted)
2464	!--------------------------------------------------------------------------------------------------!
2465
2466	RECURSIVE FUNCTION chem_emissions_locate_timestep &
2467	( this_timestamp, timestamp_array, lower_bound, upper_bound ) &
2468	RESULT ( timestamp_index )
2469
2470	!
2471	!-- Function arguments
2472	CHARACTER(LEN=*), INTENT(IN) :: this_timestamp !> target timestamp
2473	CHARACTER(LEN=512), INTENT(IN) :: timestamp_array(:) !> array of timestamps
2474
2475	INTEGER(iwp), INTENT(IN) :: lower_bound, upper_bound !> timestamp_array index bounds
2476
2477	INTEGER(iwp) :: timestamp_index !> index for most recent timestamp in timestamp_array
2478
2479	!
2480	!-- Local variables
2481	INTEGER(iwp) :: k0,km,k1 !> lower, central, and upper index bounds
2482	!
2483	!-- Assign bounds
2484	k0 = lower_bound
2485	k1 = upper_bound
2486	!
2487	!-- Make sure k1 is always not smaller than k0
2488	IF ( k0 > k1 ) THEN
2489	k0 = upper_bound
2490	k1 = lower_bound
2491	ENDIF
2492	!
2493	!-- Make sure k0 and k1 stay within array bounds by timestamp_array
2494	IF ( k0 < 1 ) k0 = 1
2495	IF ( k1 > SIZE( timestamp_array, 1 ) ) k1 = SIZE( timestamp_array, 1 )
2496	!
2497	!-- Terminate if target is contained within 2 consecutive indices otherwise calculate central bound
2498	!-- (km) and determine new index bounds for the next iteration
2499
2500	IF ( ( k1 - k0 ) > 1 ) THEN
2501	km = ( k0 + k1 ) / 2
2502	IF ( TRIM( this_timestamp ) > TRIM( timestamp_array(km) ) ) THEN
2503	k0 = km
2504	ELSE
2505	k1 = km
2506	ENDIF
2507	timestamp_index = chem_emissions_locate_timestep( this_timestamp, timestamp_array, k0, k1 )
2508	ELSE
2509	timestamp_index = k0
2510	ENDIF
2511
2512	END FUNCTION chem_emissions_locate_timestep
2513
2514
2515	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2516	!!
2517	!! END OF MODULE
2518	!!
2519	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2520
2521	END MODULE chem_emissions_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |