Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |