Home

Context Navigation

chem_emissions_mod.f90 @ 4899

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

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

Context Navigation

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

Download in other formats: