Home

Context Navigation

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

Last change on this file since 4808 was 4671, checked in by pavelkrc, 4 years ago
Radiative transfer model RTM version 4.1
Property svn:keywords set to `Id`
File size: 97.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |