module_input_chem_data.F
References to this file elsewhere.
1 !dis
2 !dis Open Source License/Disclaimer, Forecast Systems Laboratory
3 !dis NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
4 !dis
5 !dis This software is distributed under the Open Source Definition,
6 !dis which may be found at http://www.opensource.org/osd.html.
7 !dis
8 !dis In particular, redistribution and use in source and binary forms,
9 !dis with or without modification, are permitted provided that the
10 !dis following conditions are met:
11 !dis
12 !dis - Redistributions of source code must retain this notice, this
13 !dis list of conditions and the following disclaimer.
14 !dis
15 !dis - Redistributions in binary form must provide access to this
16 !dis notice, this list of conditions and the following disclaimer, and
17 !dis the underlying source code.
18 !dis
19 !dis - All modifications to this software must be clearly documented,
20 !dis and are solely the responsibility of the agent making the
21 !dis modifications.
22 !dis
23 !dis - If significant modifications or enhancements are made to this
24 !dis software, the FSL Software Policy Manager
25 !dis (softwaremgr@fsl.noaa.gov) should be notified.
26 !dis
27 !dis THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
28 !dis AND ARE FURNISHED "AS IS." THE AUTHORS, THE UNITED STATES
29 !dis GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
30 !dis AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
31 !dis OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME
32 !dis NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
33 !dis DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
34 !dis
35 !dis
36
37 !WRF:PACKAGE:IO
38
39
40 !CPP directives to control ic/bc conditions...
41 !(The directive in module_mosaic_initmixrats also needs to be set.)
42 ! CASENAME = 0 Uses Texas August 2004 case values and profiles
43 ! 1 Uses same concentrations as TX, but uses different
44 ! profiles depending on the species. (NEAQS2004 case)
45 #define CASENAME 0
46
47
48 MODULE module_input_chem_data
49
50 USE module_io_domain
51 USE module_domain
52 USE module_driver_constants
53 USE module_state_description
54 USE module_configure
55 USE module_date_time
56 USE module_wrf_error
57 USE module_timing
58 USE module_data_radm2
59 USE module_aerosols_sorgam
60 USE module_data_sorgam
61 USE module_utility
62 USE module_get_file_names
63
64
65 IMPLICIT NONE
66
67 ! REAL :: grav = 9.8104
68 REAL, PARAMETER :: mwso4 = 96.0576
69
70
71 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
72 ! Initial atmospheric chemistry profile data
73 INTEGER :: k_loop ! Used for loop index
74 INTEGER :: lo ! number of chemicals in inital profile
75 INTEGER :: logg ! number of final chemical species (nch-1)
76 INTEGER :: kx ! number of vertical levels in temp profile
77 INTEGER :: kxm1
78
79 PARAMETER( kx=16, kxm1=kx-1, logg=100, lo=34)
80
81 INTEGER, DIMENSION(logg) :: iref
82
83 REAL, DIMENSION(logg) :: fracref
84 REAL, DIMENSION(kx) :: dens
85 REAL, DIMENSION(kx+1) :: zfa
86 REAL, DIMENSION(kx+1) :: zfa_bdy
87 REAL, DIMENSION(lo ,kx) :: xl
88 REAL :: so4vaptoaer
89 DATA so4vaptoaer/.999/
90
91 CHARACTER (LEN=4), DIMENSION(logg) :: ggnam
92
93 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
94 !
95 ! The idealized profile is converted from the NALROM chemistry model.
96 !
97 ! The variable "iref" is the reference index
98 ! "fracref" is the reference fraction correpsonding to iref
99 !
100 ! For example: WRF-Chem species number 1 is SO2. iref(1) is 12, and XL(12,K)
101 ! is the profile for SO2.
102 !
103 ! Note: NALROM calculates lumped OX (XL(1) =O3+NO2+HNO3+...) and a
104 ! lumped NOX (XL(2)=NO+NO2+NO3+2N2O5+HO2NO2+HONO). But XL(33) is
105 ! strictly O3, and XL(34)=NOx=(NO+NO2 only).
106 !
107 ! Short-lived species are initialized to steady-state equilibrium -
108 ! since they are short-lived. The short-lived species within a lumped category
109 ! (Ox , NOx, or NO3+N2O5 in our case) would be renormalized to the lumped class
110 ! after the steady-state equilibrium concentrations are determined.
111 !
112 ! The following is a list of long-lived species
113 !
114 ! NAMEL( 1)='OX '
115 ! NAMEL( 2)='NOX '
116 ! NAMEL( 3)='HNO3 '
117 ! NAMEL( 4)='H2O2 '
118 ! NAMEL( 5)='CH3OOH '
119 ! NAMEL( 6)='CO '
120 ! NAMEL( 7)='ISOPRENE '
121 ! NAMEL( 8)='CH2O '
122 ! NAMEL( 9)='CH3CHO '
123 ! NAMEL(10)='PAN '
124 ! NAMEL(11)='OTHER ALKA'
125 ! NAMEL(12)='SO2 '
126 ! NAMEL(13)='BUTANE '
127 ! NAMEL(14)='ETHENE '
128 ! NAMEL(15)='PROPENE+ '
129 ! NAMEL(16)='PPN '
130 ! NAMEL(17)='MEK '
131 ! NAMEL(18)='RCHO '
132 ! NAMEL(19)='SO4= '
133 ! NAMEL(20)='MVK '
134 ! NAMEL(21)='MACR '
135 ! NAMEL(22)='HAC '
136 ! NAMEL(23)='MGLY '
137 ! NAMEL(24)='HPAN '
138 ! NAMEL(25)='MPAN '
139 ! NAMEL(26)='PROPANE '
140 ! NAMEL(27)='ACETYLENE '
141 ! NAMEL(28)='OH '
142 ! NAMEL(29)='HO2 '
143 ! NAMEL(30)='NO3+N2O5 '
144 ! NAMEL(31)='HO2NO2 '
145 ! NAMEL(32)='SUM RO2 '
146 ! NAMEL(33)='OZONE '
147 ! NAMEL(34)='NOX '
148 !
149 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
150
151
152 DATA iref/12,19,2,2,1,3,4,9,8,5,5,32,6,6,6,30,30,10,26,13,11,6,6, &
153 14,15,15,23,23,32,16,23,31,17,5*23,7,28,29,59*7/
154
155 DATA fracref/1.,1.,.75,.25,6*1.,.5,.5,6.25E-4,7.5E-4,6.25E-5,.1, &
156 .9,4*1.,8.E-3,3*1.,.5,1.,1.,.5,12*1.,59*1./
157
158 DATA ggnam/ 'SO2 ','SULF','NO2 ','NO ','O3 ','HNO3', &
159 'H2O2','ALD ','HCHO','OP1 ','OP2 ','PAA ', &
160 'ORA1','ORA2','NH3 ','N2O5','NO3 ','PAN ', &
161 'HC3 ','HC5 ','HC8 ','ETH ','CO ','OL2 ', &
162 'OLT ','OLI ','TOL ','XYL ','ACO3','TPAN', &
163 'HONO','HNO4','KET ','GLY ','MGLY','DCB ', &
164 'ONIT','CSL ','ISO ','HO ','HO2 ',59*'JUNK' /
165
166 DATA dens/ 2.738E+18, 5.220E+18, 7.427E+18, 9.202E+18, &
167 1.109E+19, 1.313E+19, 1.525E+19, 1.736E+19, &
168 1.926E+19, 2.074E+19, 2.188E+19, 2.279E+19, &
169 2.342E+19, 2.384E+19, 2.414E+19, 2.434E+19 /
170
171 ! Profile heights in meters
172 ! DATA ZFA/ 0., 85., 212., 385., 603., 960., 1430., 2010., &
173 ! 2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
174 ! 21000./
175 #if ( ! EM_CORE == 0 )
176 DATA ZFA_BDY/ 0., 85., 212., 385., 603., 960., 1430., 2010., &
177 2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
178 21000./
179
180 ! Profile heights in meters
181 DATA ZFA/ 0., 85., 212., 385., 603., 960., 1430., 2010., &
182 2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
183 21000./
184 #endif
185 #if ( ! NMM_CORE == 0 )
186
187 DATA ZFA_BDY/ 0., 85., 212., 385., 603., 960., 1430., 2010., &
188 2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
189 21000./
190
191 ! Profile pressure in hpa
192 DATA ZFA/ 100000., 98500., 98000., 96000., 94000., 90000., 85000., 75000., &
193 71000., 65000., 52000., 48000., 45000., 30000., 25000., 20000., &
194 5000./
195 #endif
196
197 !wig: To match the xl profile to the correct species, match WRF's p_<species>
198 ! flag with iref(p_<species>-1) to get the value of the first index in xl,
199 ! e.g. p_o3=6, iref(6-1)=1, so xl(1,:) is the ozone profile.
200 ! See gasprofile_init_pnnl for an explination of what height
201 ! each index represents.
202 DATA (xl(1,k_loop),k_loop=1,kx) &
203 / 1.68E-07, 1.68E-07, 5.79E-08, 5.24E-08, 5.26E-08, &
204 5.16E-08, 4.83E-08, 4.50E-08, 4.16E-08, 3.80E-08, 3.56E-08, &
205 3.35E-08, 3.15E-08, 3.08E-08, 3.06E-08, 3.00E-08/
206
207 DATA (xl(2,k_loop),k_loop=1,kx) &
208 / 4.06E-10, 4.06E-10, 2.16E-10, 1.37E-10, 9.47E-11, &
209 6.95E-11, 5.31E-11, 4.19E-11, 3.46E-11, 3.01E-11, 2.71E-11, &
210 2.50E-11, 2.35E-11, 2.26E-11, 2.20E-11, 2.16E-11/
211
212 DATA (xl(3,k_loop),k_loop=1,kx) &
213 / 9.84E-10, 9.84E-10, 5.66E-10, 4.24E-10, 3.26E-10, &
214 2.06E-10, 1.12E-10, 7.33E-11, 7.03E-11, 7.52E-11, 7.96E-11, &
215 7.56E-11, 7.27E-11, 7.07E-11, 7.00E-11, 7.00E-11/
216
217 DATA (xl(4,k_loop),k_loop=1,kx) &
218 / 8.15E-10, 8.15E-10, 8.15E-10, 8.15E-10, 8.15E-10, &
219 8.65E-10, 1.07E-09, 1.35E-09, 1.47E-09, 1.47E-09, 1.47E-09, &
220 1.47E-09, 1.45E-09, 1.43E-09, 1.40E-09, 1.38E-09/
221
222 DATA (xl(5,k_loop),k_loop=1,kx) &
223 / 4.16E-10, 4.16E-10, 4.16E-10, 4.16E-10, 4.16E-10, &
224 4.46E-10, 5.57E-10, 1.11E-09, 1.63E-09, 1.63E-09, 1.63E-09, &
225 1.63E-09, 1.61E-09, 1.59E-09, 1.57E-09, 1.54E-09/
226
227 ! CO is 70 ppbv at top, 80 throughout troposphere
228 DATA (xl(6,k_loop),k_loop=1,kx) / 7.00E-08, kxm1*8.00E-08/
229
230 DATA (xl(7,k_loop),k_loop=1,kx) &
231 / 8.33E-29, 8.33E-29, 8.33E-29, 8.33E-29, 8.33E-29, &
232 1.33E-28, 3.54E-28, 1.85E-28, 1.29E-29, 1.03E-30, 1.72E-31, &
233 7.56E-32, 1.22E-31, 2.14E-31, 2.76E-31, 2.88E-31/
234
235 DATA (xl(8,k_loop),k_loop=1,kx) &
236 / 9.17E-11, 9.17E-11, 9.17E-11, 9.17E-11, 9.17E-11, &
237 1.03E-10, 1.55E-10, 2.68E-10, 4.47E-10, 4.59E-10, 4.72E-10, &
238 4.91E-10, 5.05E-10, 5.13E-10, 5.14E-10, 5.11E-10/
239 DATA (xl(9,k_loop),k_loop=1,kx) &
240 / 7.10E-12, 7.10E-12, 7.10E-12, 7.10E-12, 7.10E-12, &
241 7.36E-12, 1.02E-11, 2.03E-11, 2.98E-11, 3.01E-11, 3.05E-11, &
242 3.08E-11, 3.08E-11, 3.06E-11, 3.03E-11, 2.99E-11/
243 DATA (xl(10,k_loop),k_loop=1,kx) &
244 / 4.00E-11, 4.00E-11, 4.00E-11, 3.27E-11, 2.51E-11, &
245 2.61E-11, 2.20E-11, 1.69E-11, 1.60E-11, 1.47E-11, 1.37E-11, &
246 1.30E-11, 1.24E-11, 1.20E-11, 1.18E-11, 1.17E-11/
247 DATA (xl(11,k_loop),k_loop=1,kx) &
248 / 1.15E-16, 1.15E-16, 2.46E-15, 2.30E-14, 1.38E-13, &
249 6.25E-13, 2.31E-12, 7.32E-12, 1.87E-11, 3.68E-11, 6.10E-11, &
250 9.05E-11, 1.22E-10, 1.50E-10, 1.70E-10, 1.85E-10/
251 DATA (xl(12,k_loop),k_loop=1,kx) &
252 / 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, &
253 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, &
254 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10/
255 DATA (xl(13,k_loop),k_loop=1,kx) &
256 / 1.26E-11, 1.26E-11, 2.02E-11, 2.50E-11, 3.02E-11, &
257 4.28E-11, 6.62E-11, 1.08E-10, 1.54E-10, 2.15E-10, 2.67E-10, &
258 3.24E-10, 3.67E-10, 3.97E-10, 4.16E-10, 4.31E-10/
259 DATA (xl(14,k_loop),k_loop=1,kx) &
260 / 1.15E-16, 1.15E-16, 2.46E-15, 2.30E-14, 1.38E-13, &
261 6.25E-13, 2.31E-12, 7.32E-12, 1.87E-11, 3.68E-11, 6.10E-11, &
262 9.05E-11, 1.22E-10, 1.50E-10, 1.70E-10, 1.85E-10/
263 DATA (xl(15,k_loop),k_loop=1,kx) &
264 / 1.00E-20, 1.00E-20, 6.18E-20, 4.18E-18, 1.23E-16, &
265 2.13E-15, 2.50E-14, 2.21E-13, 1.30E-12, 4.66E-12, 1.21E-11, &
266 2.54E-11, 4.47E-11, 6.63E-11, 8.37E-11, 9.76E-11/
267 DATA (xl(16,k_loop),k_loop=1,kx) &
268 / 1.23E-11, 1.23E-11, 1.23E-11, 1.23E-11, 1.23E-11, &
269 1.20E-11, 9.43E-12, 3.97E-12, 1.19E-12, 1.11E-12, 9.93E-13, &
270 8.66E-13, 7.78E-13, 7.26E-13, 7.04E-13, 6.88E-13/
271 DATA (xl(17,k_loop),k_loop=1,kx) &
272 / 1.43E-12, 1.43E-12, 1.43E-12, 1.43E-12, 1.43E-12, &
273 1.50E-12, 2.64E-12, 8.90E-12, 1.29E-11, 1.30E-11, 1.32E-11, &
274 1.32E-11, 1.31E-11, 1.30E-11, 1.29E-11, 1.27E-11/
275 DATA (xl(18,k_loop),k_loop=1,kx) &
276 / 3.61E-13, 3.61E-13, 3.61E-13, 3.61E-13, 3.61E-13, &
277 3.58E-13, 5.22E-13, 1.75E-12, 2.59E-12, 2.62E-12, 2.64E-12, &
278 2.66E-12, 2.65E-12, 2.62E-12, 2.60E-12, 2.57E-12/
279 DATA (xl(19,k_loop),k_loop=1,kx) &
280 / 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, &
281 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, &
282 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11/
283
284 DATA (xl(20,k_loop),k_loop=1,kx)/kx*1.E-20/
285 DATA (xl(21,k_loop),k_loop=1,kx)/kx*1.E-20/
286 DATA (xl(22,k_loop),k_loop=1,kx)/kx*1.E-20/
287 DATA (xl(23,k_loop),k_loop=1,kx)/kx*1.E-20/
288 DATA (xl(24,k_loop),k_loop=1,kx)/kx*1.E-20/
289 DATA (xl(25,k_loop),k_loop=1,kx)/kx*1.E-20/
290
291 ! Propane - Gregory PEM-West A 25 ppt median marine boundary layer
292 DATA (xl(26,k_loop),k_loop=1,kx) &
293 /5.00E-13, 1.24E-12, 2.21E-12, 3.27E-12, 4.71E-12, &
294 6.64E-12, 9.06E-12, 1.19E-11, 1.47E-11, 1.72E-11, &
295 1.93E-11, 2.11E-11, 2.24E-11, 2.34E-11, 2.42E-11, 2.48E-11/
296 ! Acetylene - Gregory PEM-West A 53 ppt median marine boundary layer
297 DATA (xl(27,k_loop),k_loop=1,kx) &
298 /1.00E-12, 2.48E-12, 4.42E-12, 6.53E-12, 9.42E-12, &
299 1.33E-11, 1.81E-11, 2.37E-11, 2.95E-11, 3.44E-11, &
300 3.85E-11, 4.22E-11, 4.49E-11, 4.69E-11, 4.84E-11, 4.95E-11/
301 ! OH
302 DATA (xl(28,k_loop),k_loop=1,kx) &
303 / 9.80E+06, 9.80E+06, 4.89E+06, 2.42E+06, 1.37E+06, &
304 9.18E+05, 7.29E+05, 6.26E+05, 5.01E+05, 4.33E+05, 4.05E+05, &
305 3.27E+05, 2.54E+05, 2.03E+05, 1.74E+05, 1.52E+05/
306 ! HO2
307 DATA (xl(29,k_loop),k_loop=1,kx) &
308 / 5.74E+07, 5.74E+07, 7.42E+07, 8.38E+07, 8.87E+07, &
309 9.76E+07, 1.15E+08, 1.34E+08, 1.46E+08, 1.44E+08, 1.40E+08, &
310 1.36E+08, 1.31E+08, 1.28E+08, 1.26E+08, 1.26E+08/
311 ! NO3+N2O5
312 DATA (xl(30,k_loop),k_loop=1,kx) &
313 / 5.52E+05, 5.52E+05, 3.04E+05, 2.68E+05, 2.32E+05, &
314 1.66E+05, 1.57E+05, 1.72E+05, 1.98E+05, 2.22E+05, 2.43E+05, &
315 2.75E+05, 3.00E+05, 3.18E+05, 3.32E+05, 3.39E+05/
316 ! HO2NO2
317 DATA (xl(31,k_loop),k_loop=1,kx) &
318 / 7.25E+07, 7.25E+07, 6.36E+07, 5.55E+07, 4.94E+07, &
319 3.66E+07, 2.01E+07, 9.57E+06, 4.75E+06, 2.37E+06, 1.62E+06, &
320 9.86E+05, 7.05E+05, 5.63E+05, 4.86E+05, 4.41E+05/
321 ! Sum of RO2 &
322 DATA (xl(32,k_loop),k_loop=1,kx) &
323 / 9.14E+06, 9.14E+06, 1.46E+07, 2.14E+07, 2.76E+07, &
324 3.62E+07, 5.47E+07, 1.19E+08, 2.05E+08, 2.25E+08, 2.39E+08, &
325 2.58E+08, 2.82E+08, 2.99E+08, 3.08E+08, 3.15E+08/
326 ! O3 <--This is not the O3 used for RADM2 or CBMZ (wig)
327 DATA (xl(33,k_loop),k_loop=1,kx) &
328 / 8.36E+11, 8.36E+11, 4.26E+11, 4.96E+11, 6.05E+11, &
329 6.93E+11, 7.40E+11, 7.74E+11, 7.82E+11, 7.75E+11, 7.69E+11, &
330 7.59E+11, 7.54E+11, 7.50E+11, 7.47E+11, 7.45E+11/
331 ! NOx (NO+NO2)
332 DATA (xl(34,k_loop),k_loop=1,kx) &
333 / 1.94E+09, 1.94E+09, 1.53E+09, 1.24E+09, 1.04E+09, &
334 8.96E+08, 7.94E+08, 7.11E+08, 6.44E+08, 6.00E+08, 5.70E+08, &
335 5.49E+08, 5.35E+08, 5.28E+08, 5.24E+08, 5.23E+08/
336
337 CONTAINS
338
339 SUBROUTINE input_ext_chem_file (si_grid)
340 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
341
342 IMPLICIT NONE
343
344 TYPE(domain) :: si_grid
345
346 INTEGER :: i , j , k, l, &
347 ids, ide, jds, jde, kds, kde, &
348 ims, ime, jms, jme, kms, kme, &
349 ips, ipe, jps, jpe, kps, kpe
350 INTEGER :: si_jday
351 INTEGER :: dat_jday,dat_year,dat_month,dat_day,dat_hour,dat_min,dat_sec
352 INTEGER :: time_loop_max , time_loop
353 INTEGER, DIMENSION(2) :: num_steps
354 INTEGER :: fid, ierr, rc
355 INTEGER :: status_next_var
356 INTEGER :: debug_level
357 INTEGER :: si_year,si_month,si_day,si_hour,si_min,si_sec
358 INTEGER :: total_time_sec , file_counter
359 #if ( ! NMM_CORE == 0 )
360 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: pint
361 REAL, ALLOCATABLE, DIMENSION(:,:) :: pdsl
362 #endif
363
364 LOGICAL :: input_from_file , need_new_file
365
366 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: si_zsigf, si_zsig
367 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ch_zsigf, ch_zsig, ozone
368 REAL :: ext_time, dat_time
369 REAL :: wgt0,wgt_time1,wgt_time2
370
371 CHARACTER (LEN=80) :: inpname, message
372 CHARACTER (LEN=19) :: date_string
373 CHARACTER (LEN=19) :: extract_date, use_date, current_date_char
374 CHARACTER*80 :: timestr
375
376 TYPE (grid_config_rec_type) :: config_flags
377 TYPE(domain) , POINTER :: null_domain, chem_grid, chgrid
378 TYPE(domain) , POINTER :: chem_grid2, chgrid2
379 ! TYPE(ESMF_Time) :: CurrTime
380
381 ! Interface block for routine that passes pointers and needs to know that they
382 ! are receiving pointers.
383
384
385 CALL model_to_grid_config_rec ( si_grid%id , model_config_rec , config_flags )
386 ! After current_date has been set, fill in the julgmt stuff.
387
388 CALL geth_julgmt ( config_flags%julyr , config_flags%julday , config_flags%gmt )
389
390 WRITE ( extract_date , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
391 model_config_rec%start_year (si_grid%id) , &
392 model_config_rec%start_month (si_grid%id) , &
393 model_config_rec%start_day (si_grid%id) , &
394 model_config_rec%start_hour (si_grid%id) , &
395 model_config_rec%start_minute(si_grid%id) , &
396 model_config_rec%start_second(si_grid%id)
397
398 write(message,'(A,A)') 'Subroutine input_chem: finding data at date/time: ',extract_date
399 CALL wrf_message ( TRIM(message) )
400
401
402 ! And here is an instance of using the information in the NAMELIST.
403
404 CALL nl_get_debug_level ( 1,debug_level )
405 CALL set_wrf_debug_level ( debug_level )
406
407
408 ! Allocated and configure the mother domain. Since we are in the nesting down
409 ! mode, we know a) we got a nest, and b) we only got 1 nest.
410
411 NULLIFY( null_domain )
412
413 CALL wrf_debug ( 100 , 'wrfchem:input_chem: calling alloc_and_configure_domain 1' )
414 ! Note that this is *not* the intended use of alloc_and_configure_domain()
415 ! It does not seem to hurt anything, yet...
416
417 ! if( si_grid%id .EQ. 1) then
418 CALL alloc_and_configure_domain ( domain_id = 1 , &
419 grid = chem_grid , &
420 parent = null_domain , &
421 kid = -1 )
422
423 ! else
424 ! CALL alloc_and_configure_domain ( domain_id = 2 , &
425 ! grid = chem_grid , &
426 ! parent = parent_grid , &
427 ! kid = 1 )
428 ! endif
429
430
431
432 CALL wrf_debug ( 100 , 'wrfchem:input_chem: set pointer for domain 1' )
433 chgrid => chem_grid
434
435 ! Get a list of available file names to try. This fills up the eligible_file_name
436 ! array with number_of_eligible_files entries. This routine issues a nonstandard
437 ! call (system).
438
439 file_counter = 1
440 need_new_file = .FALSE.
441
442 CALL unix_ls ( 'wrf_chem_input' , chem_grid%id )
443 write(message,'(A,A)') 'number of eligible files ', number_of_eligible_files
444 CALL wrf_message( TRIM(message) )
445
446 ! ! Open the input data (chemin_d01_000000) for reading.
447 ! CALL wrf_debug ( 100 , 'subroutine input_chem: calling open_r_dataset for wrfout' )
448 ! CALL construct_filename ( inpname , 'chemin' , chgrid%id, 2 , 0 , 6 )
449
450 CALL construct_filename2a (inpname , chgrid%input_chem_inname, chgrid%id , 2, extract_date)
451 write(message,'(A,A)') 'Subroutine input_chem: Opening data file ',TRIM(inpname)
452 CALL wrf_message( TRIM(message) )
453
454 CALL open_r_dataset ( fid, TRIM(inpname) , chgrid, config_flags, "DATASET=INPUT", ierr )
455
456 IF ( ierr .NE. 0 ) THEN
457 WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'subroutine chemin: error opening ',TRIM(inpname),' for reading ierr=',ierr
458 CALL WRF_ERROR_FATAL ( wrf_err_message )
459 ENDIF
460
461 ! How many data time levels in the input file?
462
463 num_steps = -1
464 time_loop_max = 0
465 CALL wrf_debug ( 100, 'subroutine input_chem: find time_loop_max' )
466
467 ! Which times are in this file, and more importantly, are any of them the
468 ! ones that we want? We need to loop over times in each files, loop
469 ! over files.
470
471 get_the_right_time : DO
472 ! CALL ext_ncd_get_next_time ( fid, date_string, status_next_var )
473 CALL wrf_get_next_time ( fid, date_string, status_next_var )
474
475 write(message,'(6A)') 'Subroutine input_chem: file date/time = ',date_string,' status = ', status_next_var
476 CALL wrf_message ( TRIM(message) )
477
478 IF ( status_next_var == 0 ) THEN
479 CALL wrf_debug ( 100 , 'input_ext_chem_file: calling close_dataset for ' // TRIM(eligible_file_name(file_counter)) )
480 CALL close_dataset ( fid , config_flags , "DATASET=INPUT" )
481 time_loop_max = time_loop_max + 1
482 IF ( time_loop_max .GT. number_of_eligible_files ) THEN
483 WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'program input_chem_data: opening too many files'
484 CALL WRF_ERROR_FATAL ( wrf_err_message )
485 END IF
486
487 IF ( time_loop_max .EQ. number_of_eligible_files ) THEN
488 num_steps(1) = time_loop_max
489 num_steps(2) = time_loop_max+1
490 use_date = date_string
491 wgt_time1 = dat_time
492
493 EXIT get_the_right_time
494 END IF
495 CYCLE get_the_right_time
496 ELSE
497 ! ELSE IF ( TRIM(date_string) .LT. TRIM(current_date_char) ) THEN
498 ! CYCLE get_the_right_time
499 ! ELSE IF ( TRIM(date_string) .EQ. TRIM(current_date_char) ) THEN
500 ! EXIT get_the_right_time
501 ! ELSE IF ( TRIM(date_string) .GT. TRIM(current_date_char) ) THEN
502 ! WRITE( wrf_err_message , FMT='(A,A,A,A,A)' ) 'Found ',TRIM(date_string),' before I found ',TRIM(current_date_char),' .'
503 ! CALL WRF_ERROR_FATAL ( wrf_err_message )
504 ! END IF
505
506 ! For now, the input date and time MUST match
507 !
508 ! Put the time check here and set num_steps
509
510 num_steps(1) = time_loop_max
511 num_steps(2) = time_loop_max+1
512 use_date = date_string
513 wgt_time1 = dat_time
514
515 EXIT get_the_right_time
516
517 ENDIF
518 if( num_steps(2) == time_loop_max ) then
519 wgt_time2 = dat_time
520 endif
521 END DO get_the_right_time
522
523 num_steps(2) = MIN(num_steps(2),time_loop_max)
524
525 ! wgt0 = (ext_time - wgt_time1) / (wgt_time2 - wgt_time1)
526 wgt0 = 0.
527
528 ! Make sure the right date and time for the chemin data has been found
529
530 write(message,'(A,A20,A,I9)') 'Subroutine input_chem: use_date, num_steps(1) = ',use_date,num_steps(1)
531 if( num_steps(1) > 0 ) then
532 write(message,'(A,A)') 'Subroutine input_chem: extracting data at date/time: ',use_date
533 CALL wrf_message ( TRIM( message ) )
534 else
535 WRITE( wrf_err_message, FMT='(A)' ) 'subroutine input_chem: error finding chemin date/time #2'
536 CALL WRF_ERROR_FATAL ( wrf_err_message )
537 endif
538
539 ! There has to be a more elegant way to get to the beginning of the file, but this will do.
540
541 CALL close_dataset ( fid , config_flags , "DATASET=INPUT" )
542
543 CALL construct_filename2a (inpname , chgrid%input_chem_inname, chgrid%id , 2, extract_date)
544 write(message,'(A,A)') 'Subroutine input_chem: Opening data file ',TRIM(inpname)
545 CALL wrf_message( TRIM(message) )
546
547 CALL open_r_dataset ( fid, TRIM(inpname) , chgrid , config_flags , "DATASET=INPUT", ierr )
548 IF ( ierr .NE. 0 ) THEN
549 WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'subroutine chemin: error opening ',TRIM(inpname),' for reading ierr=',ierr
550 CALL WRF_ERROR_FATAL ( wrf_err_message )
551 ENDIF
552
553 ! We know how many time periods to process (right now - all of them), we have the input data
554 ! (re-)opened, so we begin.
555
556 big_time_loop_thingy : DO time_loop = 1 , time_loop_max
557
558 CALL wrf_debug ( 100 , 'input_chem: calling input_history' )
559 CALL input_history ( fid , chgrid , config_flags, ierr )
560 CALL wrf_debug ( 100 , 'input_chem: back from input_history' )
561
562 IF( time_loop .EQ. num_steps(1) ) THEN
563
564 ! Get grid dimensions
565 CALL get_ijk_from_grid ( si_grid , &
566 ids, ide, jds, jde, kds, kde, &
567 ims, ime, jms, jme, kms, kme, &
568 ips, ipe, jps, jpe, kps, kpe )
569
570 ! Get scalar grid point heights
571 ALLOCATE( si_zsigf(ims:ime,kms:kme,jms:jme) )
572 ALLOCATE( ch_zsigf(ims:ime,kms:kme,jms:jme) )
573 ALLOCATE( si_zsig(ims:ime,kms:kme,jms:jme) )
574 ALLOCATE( ch_zsig(ims:ime,kms:kme,jms:jme) )
575 #if ( EM_CORE == 1 )
576 si_zsigf = (si_grid%em_ph_1 + si_grid%em_phb) / grav
577 ch_zsigf = ( chgrid%em_ph_1 + chgrid%em_phb) / grav
578 #endif
579
580 #if ( NMM_CORE == 1 )
581
582 ! Get scalar grid point heights
583 ALLOCATE( pint(ips:ipe,kps:kpe,jps:jpe) )
584 ALLOCATE( pdsl(ips:ipe,jps:jpe) )
585
586 IF(chgrid%sigma.EQ. 1)THEN
587 do j=jps,jpe
588 do i=ips,ipe
589 pdsl(i,j)=chgrid%nmm_pd(i,j)
590 ENDDO
591 ENDDO
592 ELSE
593 do j=jps,jpe
594 do i=ips,ipe
595 pdsl(i,j)=chgrid%nmm_res(i,j)*chgrid%nmm_pd(i,j)
596 enddo
597 enddO
598 ENDIF
599 !!
600 !!!!! SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
601 !
602 do j=jps,jpe
603 do k=kps,kpe
604 do i=ips,ipe
605 pint(i,k,j)=si_grid%nmm_eta1(k)*chgrid%nmm_pdtop+si_grid%nmm_eta2(k)*pdsl(i,j)+chgrid%nmm_pt
606 ch_zsigf(i,k,j)=pint(i,k,j)
607 ENDDO
608 ENDDO
609 ENDDO
610
611 CALL wrf_debug (0, message)
612
613 IF(si_grid%sigma.EQ. 1)THEN
614 do j=jps,jpe
615 do i=ips,ipe
616 pdsl(i,j)=si_grid%nmm_pd(i,j)
617 ENDDO
618 ENDDO
619 ELSE
620 do j=jps,jpe
621 do i=ips,ipe
622 pdsl(i,j)=si_grid%nmm_res(i,j)*si_grid%nmm_pd(i,j)
623 enddo
624 enddO
625 ENDIF
626 ! write(message,'(1e15.6,i6)') pdsl(ips,jps), si_grid%sigma
627 !!
628 !!!!! SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
629 !
630 do j=jps,jpe
631 do k=kps,kpe
632 do i=ips,ipe
633 pint(i,k,j)=si_grid%nmm_eta1(k)*si_grid%nmm_pdtop+si_grid%nmm_eta2(k)*pdsl(i,j)+si_grid%nmm_pt
634 ! if (i.EQ. ips .and. j .EQ. jps) then
635 ! print *,k,pint(i,k,j),si_grid%nmm_eta1(k),si_grid%nmm_pdtop,si_grid%nmm_eta2(k),pdsl(i,j),si_grid%nmm_pt
636 ! endif
637 si_zsigf(i,k,j)=pint(i,k,j)
638 ENDDO
639 ENDDO
640 ENDDO
641
642 ! write(message,'(4e15.6)') si_zsigf(1,1:4,1)
643 ! CALL wrf_error_fatal (message)
644
645 DEALLOCATE( pint); DEALLOCATE( pdsl )
646 #endif
647
648
649 do k=1,kde-1
650 si_zsig(:,k,:) = 0.5 * ( si_zsigf(:,k,:) + si_zsigf(:,k+1,:) )
651 ch_zsig(:,k,:) = 0.5 * ( ch_zsigf(:,k,:) + ch_zsigf(:,k+1,:) )
652 enddo
653 si_zsig(:,kde,:) = 0.5 * ( 3. * si_zsigf(:,kde,:) - si_zsigf(:,kde-1,:) )
654 ch_zsig(:,kde,:) = 0.5 * ( 3. * ch_zsigf(:,kde,:) - ch_zsigf(:,kde-1,:) )
655
656 ! check size of si_grid vs chgrid
657
658 IF( size(si_grid%chem,1) .NE. ime-ims+1 .OR. &
659 size(si_grid%chem,2) .NE. kme-kms+1 .OR. &
660 size(si_grid%chem,3) .NE. jme-jms+1 .OR. &
661 size(si_grid%chem,4) .NE. num_chem ) then
662
663 CALL wrf_debug (100, ' SI grid dimensions ' )
664 write(message,'(4i8.8)') size(si_grid%chem,1),size(si_grid%chem,2), &
665 size(si_grid%chem,3),size(si_grid%chem,4)
666 CALL wrf_debug (100, message)
667 CALL wrf_debug (100, ' Input data dimensions ' )
668 write(message,'(4i8.8)') ime-ims+1,kme-kms+1,jme-jms+1,num_chem
669 CALL wrf_debug (100, message)
670 write(wrf_err_message,'(A)') 'ERROR IN MODULE_INPUT_CHEM: bad dimensions in input data '
671 CALL WRF_ERROR_FATAL ( wrf_err_message )
672 ENDIF
673
674 ! Fill top level to prevent spurious interpolation results (no extrapolation)
675 chgrid%chem(:,kde,:,:) = chgrid%chem(:,kde-1,:,:)
676
677 ! Interpolate the chemistry data to the SI grid (holding place for time interpolation)
678
679 call vinterp_chem(ims, ime, jms, jme, kms, kme, kme, num_chem, ch_zsig, si_zsig, &
680 chgrid%chem, si_grid%chem, .false.)
681
682 if(wgt0 == 0.) EXIT big_time_loop_thingy
683 ENDIF
684
685 IF( time_loop .EQ. num_steps(2) ) THEN
686
687 ! ! input chemistry sigma levels
688 ! ch_zsigf = ( chgrid%em_ph_1 + chgrid%em_phb) / grav
689 ! do k=1,kde-1
690 ! ch_zsig(:,k,:) = 0.5 * ( ch_zsigf(:,k,:) + ch_zsigf(:,k+1,:) )
691 ! enddo
692 ! ch_zsig(:,kde,:) = 0.5 * ( 3. * ch_zsigf(:,kde,:) - ch_zsigf(:,kde-1,:) )
693
694 ! ! Fill top level to prevent spurious interpolation results (no extrapolation)
695 ! chgrid%chem(:,kde,:,:) = chgrid%chem(:,kde-1,:,:)
696
697 ! ! Interpolate the chemistry data to the temp chgrid2 structure
698
699 ! call vinterp_chem(ims, ime, jms, jme, kms, kme, kme, num_chem, ch_zsig, si_zsig, &
700 ! chgrid%chem, chgrid2%chem, .false.)
701
702 ! ! use linear interpolation in time to get new chem arrays
703 ! si_grid%chem = (1. - wgt0) * si_grid%chem + &
704 ! wgt0 * chgrid2%chem
705
706 DEALLOCATE( si_zsigf); DEALLOCATE( si_zsig )
707 DEALLOCATE( ch_zsigf); DEALLOCATE( ch_zsig )
708
709 EXIT big_time_loop_thingy
710 ENDIF
711 END DO big_time_loop_thingy
712
713 ! Check for errors in chemin data set
714
715 do l=2,num_chem
716 do j=jds,jde
717 do k=kds,kde
718 do i=ids,ide
719 si_grid%chem(i,k,j,l) = MAX(si_grid%chem(i,k,j,l),epsilc)
720 enddo
721 enddo
722 enddo
723 enddo
724
725
726 ! Close chemin data set
727 CALL close_dataset ( fid , config_flags , "DATASET=INPUT" )
728
729 ! free memory
730 CALL domain_destroy( chem_grid )
731
732 CALL wrf_debug ( 100,' input_ext_chem_data: exit subroutine ')
733
734 RETURN
735
736 END SUBROUTINE input_ext_chem_file
737 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
738 SUBROUTINE vinterp_chem(nx1, nx2, ny1, ny2, nz1, nz_in, nz_out, nch, z_in, z_out, &
739 data_in, data_out, extrapolate)
740
741 ! Interpolates columns of chemistry data from one set of height surfaces to
742 ! another
743
744 INTEGER, INTENT(IN) :: nx1, nx2
745 INTEGER, INTENT(IN) :: ny1, ny2
746 INTEGER, INTENT(IN) :: nz1
747 INTEGER, INTENT(IN) :: nz_in
748 INTEGER, INTENT(IN) :: nz_out
749 INTEGER, INTENT(IN) :: nch
750 REAL, INTENT(IN) :: z_in (nx1:nx2,nz1:nz_in ,ny1:ny2)
751 REAL, INTENT(IN) :: z_out(nx1:nx2,nz1:nz_out,ny1:ny2)
752 REAL, INTENT(IN) :: data_in (nx1:nx2,nz1:nz_in ,ny1:ny2,nch)
753 REAL, INTENT(OUT) :: data_out(nx1:nx2,nz1:nz_out,ny1:ny2,nch)
754 LOGICAL, INTENT(IN) :: extrapolate
755
756 INTEGER :: i,j,l
757 INTEGER :: k,kk
758 REAL :: desired_z
759 REAL :: dvaldz
760 REAL :: wgt0
761
762 ! Loop over the number of chemical species
763 chem_loop: DO l = 2, nch
764
765 data_out(:,:,:,l) = -99999.9
766
767 DO j = ny1, ny2
768 DO i = nx1, nx2
769
770 output_loop: DO k = nz1, nz_out
771 #if ( EM_CORE == 1 )
772
773 desired_z = z_out(i,k,j)
774 IF (desired_z .LT. z_in(i,1,j)) THEN
775
776 IF ((desired_z - z_in(i,1,j)).LT. 0.0001) THEN
777 data_out(i,k,j,l) = data_in(i,1,j,l)
778 ELSE
779 IF (extrapolate) THEN
780 ! Extrapolate downward because desired height level is below
781 ! the lowest level in our input data. Extrapolate using simple
782 ! 1st derivative of value with respect to height for the bottom 2
783 ! input layers.
784
785 ! Add a check to make sure we are not using the gradient of
786 ! a very thin layer
787
788 IF ( (z_in(i,1,j) - z_in(i,2,j)) .GT. 0.001) THEN
789 dvaldz = (data_in(i,1,j,l) - data_in(i,2,j,l)) / &
790 (z_in(i,1,j) - z_in(i,2,j) )
791 ELSE
792 dvaldz = (data_in(i,1,j,l) - data_in(i,3,j,l)) / &
793 (z_in(i,1,j) - z_in(i,3,j) )
794 ENDIF
795 data_out(i,k,j,l) = MAX( data_in(i,1,j,l) + &
796 dvaldz * (desired_z-z_in(i,1,j)), 0.)
797 ELSE
798 data_out(i,k,j,l) = data_in(i,1,j,l)
799 ENDIF
800 ENDIF
801 ELSE IF (desired_z .GT. z_in(i,nz_in,j)) THEN
802 IF ( (z_in(i,nz_in,j) - desired_z) .LT. 0.0001) THEN
803 data_out(i,k,j,l) = data_in(i,nz_in,j,l)
804 ELSE
805 IF (extrapolate) THEN
806 ! Extrapolate upward
807 IF ( (z_in(i,nz_in-1,j)-z_in(i,nz_in,j)) .GT. 0.0005) THEN
808 dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-1,j,l)) / &
809 (z_in(i,nz_in,j) - z_in(i,nz_in-1,j))
810 ELSE
811 dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-2,j,l)) / &
812 (z_in(i,nz_in,j) - z_in(i,nz_in-2,j))
813 ENDIF
814 data_out(i,k,j,l) = MAX( data_in(i,nz_in,j,l) + &
815 dvaldz * (desired_z-z_in(i,nz_in,j)), 0.)
816 ELSE
817 data_out(i,k,j,l) = data_in(i,nz_in,j,l)
818 ENDIF
819 ENDIF
820 ELSE
821 ! We can trap between two levels and linearly interpolate
822
823 input_loop: DO kk = 1, nz_in-1
824 IF (desired_z .EQ. z_in(i,kk,j) )THEN
825 data_out(i,k,j,l) = data_in(i,kk,j,l)
826 EXIT input_loop
827 ELSE IF ( (desired_z .GT. z_in(i,kk,j)) .AND. &
828 (desired_z .LT. z_in(i,kk+1,j)) ) THEN
829 wgt0 = (desired_z - z_in(i,kk+1,j)) / &
830 (z_in(i,kk,j)-z_in(i,kk+1,j))
831 data_out(i,k,j,l) = MAX( wgt0*data_in(i,kk,j,l) + &
832 (1.-wgt0)*data_in(i,kk+1,j,l), 0.)
833 EXIT input_loop
834 ENDIF
835 ENDDO input_loop
836
837 ENDIF
838 #endif
839 #if ( NMM_CORE == 1 )
840
841 desired_z = z_out(i,k,j)
842 IF (desired_z .GT. z_in(i,1,j)) THEN
843
844 IF ((desired_z - z_in(i,1,j)).GT. 0.0001) THEN
845 data_out(i,k,j,l) = data_in(i,1,j,l)
846 ELSE
847 IF (extrapolate) THEN
848 ! Extrapolate upward because desired pressure level is above
849 ! the highest level in our input data. Extrapolate using simple
850 ! 1st derivative of value with respect to height for the bottom 2
851 ! input layers.
852
853 ! Add a check to make sure we are not using the gradient of
854 ! a very thin layer
855
856 IF ( (z_in(i,1,j) - z_in(i,2,j)) .LT. 0.001) THEN
857 dvaldz = (data_in(i,2,j,l) - data_in(i,1,j,l)) / &
858 (z_in(i,2,j) - z_in(i,1,j) )
859 ELSE
860 dvaldz = (data_in(i,3,j,l) - data_in(i,1,j,l)) / &
861 (z_in(i,3,j) - z_in(i,1,j) )
862 ENDIF
863 data_out(i,k,j,l) = MAX( data_in(i,1,j,l) + &
864 dvaldz * (desired_z-z_in(i,1,j)), 0.)
865 ELSE
866 data_out(i,k,j,l) = data_in(i,1,j,l)
867 ENDIF
868 ENDIF
869 ELSE IF (desired_z .LT. z_in(i,nz_in,j)) THEN
870 IF ( (z_in(i,nz_in,j) - desired_z) .LT. 0.0001) THEN
871 data_out(i,k,j,l) = data_in(i,nz_in,j,l)
872 ELSE
873 IF (extrapolate) THEN
874 ! Extrapolate upward
875 IF ( (z_in(i,nz_in-1,j)-z_in(i,nz_in,j)) .LT. 0.0005) THEN
876 dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-1,j,l)) / &
877 (z_in(i,nz_in,j) - z_in(i,nz_in-1,j))
878 ELSE
879 dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-2,j,l)) / &
880 (z_in(i,nz_in,j) - z_in(i,nz_in-2,j))
881 ENDIF
882 data_out(i,k,j,l) = MAX( data_in(i,nz_in,j,l) + &
883 dvaldz * (z_in(i,nz_in,j) - desired_z), 0.)
884 ELSE
885 data_out(i,k,j,l) = data_in(i,nz_in,j,l)
886 ENDIF
887 ENDIF
888 ELSE
889 ! We can trap between two levels and linearly interpolate
890
891 input_loop: DO kk = 1, nz_in-1
892 IF (desired_z .EQ. z_in(i,kk,j) )THEN
893 data_out(i,k,j,l) = data_in(i,kk,j,l)
894 EXIT input_loop
895 ELSE IF ( (desired_z .LT. z_in(i,kk,j)) .AND. &
896 (desired_z .GT. z_in(i,kk+1,j)) ) THEN
897 wgt0 = (desired_z - z_in(i,kk+1,j)) / &
898 (z_in(i,kk,j)-z_in(i,kk+1,j))
899 data_out(i,k,j,l) = MAX( wgt0*data_in(i,kk,j,l) + &
900 (1.-wgt0)*data_in(i,kk+1,j,l), 0.)
901 EXIT input_loop
902 ENDIF
903 ENDDO input_loop
904
905 ENDIF
906 #endif
907 ENDDO output_loop
908 ENDDO
909 ENDDO
910 ENDDO chem_loop
911
912 RETURN
913 END SUBROUTINE vinterp_chem
914 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
915 SUBROUTINE input_chem_profile (si_grid)
916
917 IMPLICIT NONE
918
919 TYPE(domain) :: si_grid
920
921 INTEGER :: i , j , k, &
922 ids, ide, jds, jde, kds, kde, &
923 ims, ime, jms, jme, kms, kme, &
924 ips, ipe, jps, jpe, kps, kpe
925 INTEGER :: fid, ierr, numgas
926 INTEGER :: debug_level
927
928 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: si_zsigf, si_zsig
929
930 #if ( ! NMM_CORE == 0 )
931 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: pint
932 REAL, ALLOCATABLE, DIMENSION(:,:) :: pdsl
933 #endif
934
935 CHARACTER (LEN=80) :: inpname, message
936
937 write(message,'(A)') 'Subroutine input_chem_profile: '
938 CALL wrf_message ( TRIM(message) )
939
940 ! And here is an instance of using the information in the NAMELIST.
941
942 CALL nl_get_debug_level ( 1,debug_level )
943 CALL set_wrf_debug_level ( debug_level )
944
945 ! Get grid dimensions
946 CALL get_ijk_from_grid ( si_grid , &
947 ids, ide, jds, jde, kds, kde, &
948 ims, ime, jms, jme, kms, kme, &
949 ips, ipe, jps, jpe, kps, kpe )
950
951 ! Get scalar grid point heights
952 ALLOCATE( si_zsigf(ims:ime,kms:kme,jms:jme) )
953 ALLOCATE( si_zsig(ims:ime,kms:kme,jms:jme) )
954
955 #if ( ! EM_CORE == 0 )
956 write(message,'(A)') 'WRF_EM_CORE '
957 si_zsigf = (si_grid%em_ph_1 + si_grid%em_phb) / grav
958 #endif
959 #if ( ! NMM_CORE == 0 )
960 ! Get scalar grid point heights
961 ALLOCATE( pint(ims:ime,kms:kme,jms:jme) )
962 ALLOCATE( pdsl(ims:ime,jms:jme) )
963
964 write(message,'(A)') 'WRF_NMM_CORE '
965 CALL wrf_message ( message )
966
967 IF(si_grid%sigma.EQ. 1)THEN
968 do j=jps,jpe
969 do i=ips,ipe
970 pdsl(i,j)=si_grid%nmm_pd(i,j)
971 ENDDO
972 ENDDO
973 ELSE
974 do j=jps,jpe
975 do i=ips,ipe
976 pdsl(i,j)=si_grid%nmm_res(i,j)*si_grid%nmm_pd(i,j)
977 enddo
978 enddO
979 ENDIF
980 !!
981 !!***
982 !!
983 !!
984 !!!!! SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
985 !
986 ! print *,' ips=',ips,' ipe=',ipe
987 ! print *,' jps=',jps,' jpe=',jpe
988 ! print *,' kps=',kps,' kpe=',kpe
989 ! print *,' sigma=',si_grid%sigma
990 ! print *,' pdtop=',si_grid%nmm_pdtop,' pt=',si_grid%nmm_pt
991
992 do j=jps,jpe
993 do k=kps,kpe
994 do i=ips,ipe
995 pint(i,k,j)=si_grid%nmm_eta1(k)*si_grid%nmm_pdtop+si_grid%nmm_eta2(k)*pdsl(i,j)+si_grid%nmm_pt
996 si_zsigf(i,k,j)=pint(i,k,j)
997 ENDDO
998 ENDDO
999 ENDDO
1000 ! do k=kps,kpe
1001 ! print *,k,pint(1,k,1),si_grid%nmm_eta1(k),si_grid%nmm_pdtop,si_grid%nmm_eta2(k),pdsl(1,1),si_grid%nmm_pt
1002 ! enddo
1003 !
1004 ! si_zsigf = si_grid%z
1005 #endif
1006
1007 ! si_zsigf = (si_grid%em_ph_1 + si_grid%em_phb) / grav
1008
1009 do k=1,kde-1
1010 si_zsig(:,k,:) = 0.5 * ( si_zsigf(:,k,:) + si_zsigf(:,k+1,:) )
1011 enddo
1012 si_zsig(:,kde,:) = 0.5 * ( 3. * si_zsigf(:,kde,:) - si_zsigf(:,kde-1,:) )
1013
1014 ! An alternative ozone profile option for initialization
1015 if( si_grid%gas_ic_opt == GAS_IC_PNNL ) &
1016 call gasprofile_init_pnnl
1017
1018 ! Determine the index of the last gas species
1019 numgas = get_last_gas(si_grid%chem_opt)
1020
1021 ! Interpolate the chemistry data to the SI grid. These values should typically
1022 ! be set to match the values in bdy_chem_value_tracer so that the boundaries
1023 ! and interior match each other.
1024 IF ( si_grid%chem_opt == CHEM_TRACER ) THEN
1025 si_grid%chem(ims:ime,kms:kme,jms:jme,1:numgas) = 0.0001
1026 ! si_grid%chem(ims:ime,kms:kme,jms:jme,p_so2) = 0.0001
1027 si_grid%chem(ims:ime,kms:kme,jms:jme,p_co ) = 0.08
1028 ELSE
1029 CALL make_chem_profile (ims, ime, jms, jme, kms, kme, num_chem, numgas, &
1030 si_zsig, si_grid%chem)
1031 END IF
1032
1033 CALL wrf_debug ( 100,' input_chem_profile: exit subroutine ')
1034
1035 DEALLOCATE( si_zsigf ); DEALLOCATE( si_zsig )
1036 #if ( ! NMM_CORE == 0 )
1037 DEALLOCATE( pdsl ); DEALLOCATE( pint )
1038 #endif
1039 RETURN
1040
1041 END SUBROUTINE input_chem_profile
1042 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1043 SUBROUTINE make_chem_profile ( nx1, nx2, ny1, ny2, nz1, nz2, nch, numgas, &
1044 zgrid, chem )
1045 IMPLICIT NONE
1046
1047 INTEGER, INTENT(IN) :: nx1, ny1, nz1
1048 INTEGER, INTENT(IN) :: nx2, ny2, nz2
1049 INTEGER, INTENT(IN) :: nch, numgas
1050 !REAL, INTENT(IN), DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2) :: zgrid
1051 REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2) :: zgrid
1052
1053 CHARACTER (LEN=80) :: message
1054 INTEGER :: i, j, k, l, is
1055
1056 REAL, DIMENSION(nx1:nx2,nz1:kx ,ny1:ny2,lo+1):: chprof
1057 REAL, DIMENSION(nx1:nx2,nz1:kx ,ny1:ny2) :: zprof
1058
1059 REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2,nch) :: chem
1060 REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2,lo ) :: stor
1061 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1062 !
1063 ! Check the number of species
1064
1065 if( nch .NE. num_chem) then
1066 message = ' Input_chem_profile: wrong number of chemical species'
1067 ! CALL WRF_ERROR_FATAL ( message )
1068 endif
1069
1070 ! Vertically flip the chemistry data as it is given top down and
1071 ! heights are bottom up. Fill temp 3D chemical and profile array,
1072 ! keep chem slot 1 open as vinterp_chem assumes there is no data.
1073 DO j=ny1,ny2
1074 DO k= 1,kx
1075 DO i=nx1,nx2
1076 chprof(i,k,j,2:lo+1) = xl(1:lo,kx-k+1)
1077 zprof(i,k,j) = 0.5*(zfa(k)+zfa(k+1))
1078 ENDDO
1079 ENDDO
1080 ENDDO
1081 !
1082 ! return xl to previous value for next time...
1083 ! 34 chemicals (lo), 16 vertical levels (kx)
1084 ! DO i=lo-6,lo
1085 ! xl(i,1:kx)=xl(i,1:kx)*dens(1:kx)
1086 ! ENDDO
1087
1088 ! Change number concentrations to mixing ratios for short-lived NALROM species
1089 do k=1,kx
1090 chprof(:,k,:,lo-5:lo+1) = chprof(:,k,:,lo-5:lo+1)/dens(k)
1091 end do
1092
1093 ! Interpolate temp 3D chemical and profile array to WRF grid
1094 call vinterp_chem(nx1, nx2, ny1, ny2, nz1, kx, nz2, lo, zprof, zgrid, &
1095 chprof, chem, .false.)
1096
1097 ! place interpolated data into temp storage array
1098 stor(nx1:nx2,nz1:nz2,ny1:ny2,1:lo) = chem(nx1:nx2,nz1:nz2,ny1:ny2,2:lo+1)
1099
1100 ! Here is where the chemistry profile is constructed
1101 !chem(:,:,:,1) = stor(:,:,:,1) * 0.
1102 chem(nx1:nx2,nz1:nz2,ny1:ny2,1) = -999.
1103
1104 ! DO l=2,nch
1105 DO l=2,numgas
1106 is=iref(l-1)
1107 DO j=ny1,ny2
1108 DO k=nz1,nz2
1109 DO i=nx1,nx2
1110 chem(i,k,j,l)=fracref(l-1)*stor(i,k,j,is)*1.E6
1111 ENDDO
1112 ENDDO
1113 ENDDO
1114 ENDDO
1115
1116 RETURN
1117 END SUBROUTINE make_chem_profile
1118 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1119 !
1120 ! this is a kludge routine as of now....
1121 !
1122 SUBROUTINE bdy_chem_value_sorgam (chem, z, nch, config_flags, &
1123 alt,convfac,g)
1124 USE module_data_sorgam
1125 USE module_aerosols_sorgam
1126 IMPLICIT NONE
1127
1128 REAL, intent(OUT) :: chem
1129 REAL, intent(IN) :: z ! 3D height array
1130 INTEGER, intent(IN) :: nch ! index number of chemical species
1131 REAL, INTENT(IN ) :: alt, convfac
1132 real, INTENT (IN) :: g
1133 TYPE (grid_config_rec_type), intent(in) :: config_flags
1134
1135 INTEGER :: i, k, l
1136 REAL, DIMENSION(lo+1,1:kx):: cprof ! chemical profile, diff. index order
1137
1138 REAL, DIMENSION(1:kx):: zprof
1139 REAL, DIMENSION(lo ) :: stor
1140 REAL :: wgt0
1141
1142 real :: chemsulf_radm,chem_so4aj,chem_so4ai
1143 real tempfac
1144 REAL :: splitfac
1145 !between gas and aerosol phase
1146
1147 !factor for splitting initial conc. of SO4
1148 !3rd moment i-mode [3rd moment/m^3]
1149 REAL :: m3nuc
1150 !3rd MOMENT j-mode [3rd moment/m^3]
1151 REAL :: m3acc
1152 ! REAL ESN36
1153 REAL :: m3cor
1154 DATA splitfac/.98/
1155
1156 !
1157 ! method for bc calculation is determined by aer_bc_opt
1158 !
1159 if (config_flags%aer_bc_opt == AER_BC_PNNL) then
1160 call sorgam_set_aer_bc_pnnl( chem, z, nch )
1161 return
1162 else if (config_flags%aer_bc_opt == AER_BC_DEFAULT) then
1163 continue
1164 else
1165 call wrf_error_fatal( &
1166 "bdy_chem_value_sorgam -- unable to parse aer_bc_opt" )
1167 end if
1168
1169 ! do default calculation of sorgam aerosol bc values
1170 chem=conmin
1171 ! tempfac=(t+t0)*((p+pb)/p1000mb)**rcp
1172 ! convfac=(p+pb)/rgasuniv/tempfac
1173 !
1174 !--- units for advection....
1175 !
1176 if(nch.eq.p_nu0)chem=1.e8*alt
1177 if(nch.eq.p_ac0)chem=1.e8*alt
1178 if(nch.eq.p_nh4aj)chem=10.e-5*alt
1179 if(nch.eq.p_nh4ai)chem=10.e-5*alt
1180 if(nch.eq.p_no3aj)chem=10.e-5*alt
1181 if(nch.eq.p_no3ai)chem=10.e-5*alt
1182 !
1183 ! recalculate sulf profile for aerosols
1184 !
1185 if ( nch .eq. p_so4aj.or.nch.eq.p_so4ai &
1186 .or.nch .eq. p_nu0 .or.nch.eq.p_ac0 &
1187 .or.nch .eq. p_corn ) then
1188
1189 ! Vertically flip the chemistry data as it is given top down
1190 ! and heights in zfa are bottom up
1191 ! Fill chemical profile array cprof
1192 ! Keep chem slot 1 open as vinterp_chem assumes there is no data
1193 ! (this isn't really needed in this subr)
1194 ! Convert species 28-34 (lo-6:lo) from (molecules/cm3) to (mol/mol)
1195 DO k = 1,kx
1196 zprof(k) = 0.5*(zfa_bdy(k)+zfa_bdy(k+1))
1197 DO l = 1,lo-7
1198 cprof(l+1,k) = xl(l,kx+1-k)
1199 END DO
1200 ! Fix number concentrations to mixing ratios for short-lived NALROM species
1201 DO l = lo-6,lo
1202 cprof(l+1,k) = xl(l,kx+1-k)/dens(kx+1-k)
1203 ENDDO
1204 ENDDO
1205
1206 ! Interpolate temp 1D chemical profile array to WRF field
1207 IF (z .LT. zprof(1)) THEN
1208 stor(1:lo) = cprof(2:lo+1,1)
1209 ELSE IF (z .GT. zprof(kx)) THEN
1210 stor(1:lo) = cprof(2:lo+1,kx)
1211 ELSE
1212 ! We can trap between two levels and linearly interpolate
1213 input_loop: DO k = 1, kx-1
1214 IF (z .EQ. zprof(k) )THEN
1215 stor(1:lo) = cprof(2:lo+1,k)
1216 EXIT input_loop
1217 ELSE IF ( (z .GT. zprof(k)) .AND. &
1218 (z .LT. zprof(k+1)) ) THEN
1219 wgt0 = (z - zprof(k+1)) / &
1220 (zprof(k) - zprof(k+1))
1221 stor(1:lo) = MAX( wgt0 *cprof(2:lo+1,k ) + &
1222 (1.-wgt0)*cprof(2:lo+1,k+1), 0.)
1223 EXIT input_loop
1224 ENDIF
1225 ENDDO input_loop
1226 ENDIF
1227
1228 ! Here is where the chemistry value is constructed
1229 chemsulf_radm = fracref(p_sulf-1)*stor( iref(p_sulf-1) )*1.E6
1230 !
1231 ! now have sulf
1232 !
1233 chem_so4aj=chemsulf_radm*CONVFAC*MWSO4*splitfac*so4vaptoaer
1234 chem_so4ai=chemsulf_radm*CONVFAC*MWSO4*(1.-splitfac)*so4vaptoaer
1235 if(nch.eq.p_so4aj)chem=chem_so4aj*alt
1236 if(nch.eq.p_so4ai)chem=chem_so4ai*alt
1237 m3nuc=so4fac*chem_so4ai+conmin*(nh4fac+no3fac+orgfac*9+2*anthfac)
1238 m3acc=so4fac*chem_so4aj+conmin*(nh4fac+no3fac+orgfac*9+2*anthfac)
1239 m3cor=conmin*(soilfac+seasfac+anthfac)
1240 !
1241 ! compute values for aerosol input data
1242 !
1243 if(nch.eq.p_nu0.or.nch.eq.p_ac0.or.nch.eq.p_corn)then
1244 xxlsgn = log(sginin)
1245 xxlsga = log(sginia)
1246 xxlsgc = log(sginic)
1247
1248 l2sginin = xxlsgn**2
1249 l2sginia = xxlsga**2
1250 l2sginic = xxlsgc**2
1251
1252 en1 = exp(0.125*l2sginin)
1253 ea1 = exp(0.125*l2sginia)
1254 ec1 = exp(0.125*l2sginic)
1255
1256
1257 esn04 = en1**4
1258 esa04 = ea1**4
1259 esc04 = ec1**4
1260
1261 esn05 = esn04*en1
1262 esa05 = esa04*ea1
1263
1264 esn08 = esn04*esn04
1265 esa08 = esa04*esa04
1266 esc08 = esc04*esc04
1267
1268 esn09 = esn04*esn05
1269 esa09 = esa04*esa05
1270
1271 esn12 = esn04*esn04*esn04
1272 esa12 = esa04*esa04*esa04
1273 esc12 = esc04*esc04*esc04
1274
1275 esn16 = esn08*esn08
1276 esa16 = esa08*esa08
1277 esc16 = esc08*esc08
1278
1279 esn20 = esn16*esn04
1280 esa20 = esa16*esa04
1281 esc20 = esc16*esc04
1282
1283 esn24 = esn12*esn12
1284 esa24 = esa12*esa12
1285 esc24 = esc12*esc12
1286
1287 esn25 = esn16*esn09
1288 esa25 = esa16*esa09
1289
1290 esn28 = esn20*esn08
1291 esa28 = esa20*esa08
1292 esc28 = esc20*esc08
1293
1294
1295 esn32 = esn16*esn16
1296 esa32 = esa16*esa16
1297 esc32 = esc16*esc16
1298
1299 esn36 = esn16*esn20
1300 esa36 = esa16*esa20
1301 esc36 = esc16*esc20
1302 endif
1303 !
1304 ! Units are something like number concentration
1305 !
1306 if(nch.eq.p_nu0)chem=m3nuc/((dginin**3)*esn36)*alt
1307 if(nch.eq.p_ac0)chem=m3acc/((dginia**3)*esa36)*alt
1308 if(nch.eq.p_corn)chem=m3cor/((dginic**3)*esc36)*alt
1309 endif
1310
1311
1312 END SUBROUTINE bdy_chem_value_sorgam
1313
1314 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1315 SUBROUTINE bdy_chem_value_tracer ( chem, nch )
1316
1317 ! This subroutine is called to set the boundary values of chemistry
1318 ! species when chem_opt==CHEM_TRACER. Typically, the boundary values
1319 ! here should be set to match those in input_chem_profile so that the
1320 ! interior and boundary values are the same.
1321 ! William.Gustafson@pnl.gov; 16-Jun-2005
1322
1323 IMPLICIT NONE
1324
1325 REAL, intent(OUT) :: chem
1326 INTEGER, intent(IN) :: nch ! index number of chemical species
1327 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1328
1329 if( nch .ne. p_co ) then
1330 chem = 0.0001
1331 else if( nch == p_co ) then
1332 chem = 0.08
1333 else
1334 chem = conmin
1335 end if
1336
1337 END SUBROUTINE bdy_chem_value_tracer
1338 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1339
1340 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1341 SUBROUTINE bdy_chem_value_racm ( chem, z, nch, numgas,p_co2 )
1342
1343 IMPLICIT NONE
1344
1345 REAL, intent(OUT) :: chem
1346 REAL, intent(IN) :: z ! 3D height array
1347 INTEGER, intent(IN) :: nch,p_co2 ! index number of chemical species
1348 INTEGER, intent(IN) :: numgas ! index number of last gas species
1349
1350 INTEGER :: i, k, irefcur
1351
1352 REAL, DIMENSION(kx):: cprof ! chemical profile, diff. index order
1353
1354 REAL, DIMENSION(1:kx):: zprof
1355 REAL :: stor
1356 REAL :: wgt0
1357
1358 CHARACTER (LEN=80) :: message
1359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1360
1361 ! Check the number of species
1362 ! if((nch-1).gt.logg)return
1363 if (nch.eq.p_co2)then
1364 chem=370.
1365 return
1366 endif
1367 if (nch.eq.p_co2+1)then
1368 chem=1.7
1369 return
1370 endif
1371 if (nch.ge.p_co2+2)return
1372
1373
1374 ! if( nch .GT. logg+1) then
1375 if( nch .GT. numgas) then
1376 message = ' Input_chem_profile: wrong number of chemical species'
1377 return
1378 ! CALL WRF_ERROR_FATAL ( message )
1379 endif
1380
1381 ! Vertically flip the chemistry data as it is given top down
1382 ! and heights in zfa are bottom up
1383 ! Fill 1D chemical profile array cprof
1384 ! Convert species 28-34 (lo-6:lo) from (molecules/cm3) to (mol/mol)
1385 irefcur = iref(nch-1)
1386 DO k = 1,kx
1387 zprof(k) = 0.5*(zfa_bdy(k)+zfa_bdy(k+1))
1388 if (irefcur .lt. lo-6) then
1389 cprof(k) = xl(irefcur,kx+1-k)
1390 else
1391 cprof(k) = xl(irefcur,kx+1-k)/dens(kx+1-k)
1392 end if
1393 ENDDO
1394
1395 ! Interpolate temp 3D chemical profile array to WRF field
1396 IF (z .LT. zprof(1)) THEN
1397 stor = cprof(1)
1398 ELSE IF (z .GT. zprof(kx)) THEN
1399 stor = cprof(kx)
1400 ELSE
1401 ! We can trap between two levels and linearly interpolate
1402 input_loop: DO k = 1, kx-1
1403 IF (z .EQ. zprof(k) )THEN
1404 stor = cprof(k)
1405 EXIT input_loop
1406 ELSE IF ( (z .GT. zprof(k)) .AND. &
1407 (z .LT. zprof(k+1)) ) THEN
1408 wgt0 = (z - zprof(k+1)) / &
1409 (zprof(k) - zprof(k+1))
1410 stor = MAX( wgt0 *cprof(k ) + &
1411 (1.-wgt0)*cprof(k+1), 0.)
1412 EXIT input_loop
1413 ENDIF
1414 ENDDO input_loop
1415 ENDIF
1416
1417 ! Here is where the chemistry value is constructed
1418 chem = fracref(nch-1)*stor*1.E6
1419
1420 ! special code for sulfate/h2so4
1421 if(nch.eq.p_sulf.and.p_nu0.gt.1)then
1422 chem=chem*(1.-so4vaptoaer)
1423 endif
1424
1425 RETURN
1426 END SUBROUTINE bdy_chem_value_racm
1427 SUBROUTINE bdy_chem_value ( chem, z, nch, numgas )
1428
1429 IMPLICIT NONE
1430
1431 REAL, intent(OUT) :: chem
1432 REAL, intent(IN) :: z ! 3D height array
1433 INTEGER, intent(IN) :: nch ! index number of chemical species
1434 INTEGER, intent(IN) :: numgas ! index number of last gas species
1435
1436 INTEGER :: i, k, irefcur
1437
1438 REAL, DIMENSION(kx):: cprof ! chemical profile, diff. index order
1439
1440 REAL, DIMENSION(1:kx):: zprof
1441 REAL :: stor
1442 REAL :: wgt0
1443
1444 CHARACTER (LEN=80) :: message
1445 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1446
1447 ! Check the number of species
1448 ! if((nch-1).gt.logg)return
1449
1450 ! if( nch .GT. logg+1) then
1451 if( nch .GT. numgas) then
1452 message = ' Input_chem_profile: wrong number of chemical species'
1453 return
1454 ! CALL WRF_ERROR_FATAL ( message )
1455 endif
1456
1457 ! Vertically flip the chemistry data as it is given top down
1458 ! and heights in zfa are bottom up
1459 ! Fill 1D chemical profile array cprof
1460 ! Convert species 28-34 (lo-6:lo) from (molecules/cm3) to (mol/mol)
1461 irefcur = iref(nch-1)
1462 DO k = 1,kx
1463 zprof(k) = 0.5*(zfa_bdy(k)+zfa_bdy(k+1))
1464 if (irefcur .lt. lo-6) then
1465 cprof(k) = xl(irefcur,kx+1-k)
1466 else
1467 cprof(k) = xl(irefcur,kx+1-k)/dens(kx+1-k)
1468 end if
1469 ENDDO
1470
1471 ! Interpolate temp 3D chemical profile array to WRF field
1472 IF (z .LT. zprof(1)) THEN
1473 stor = cprof(1)
1474 ELSE IF (z .GT. zprof(kx)) THEN
1475 stor = cprof(kx)
1476 ELSE
1477 ! We can trap between two levels and linearly interpolate
1478 input_loop: DO k = 1, kx-1
1479 IF (z .EQ. zprof(k) )THEN
1480 stor = cprof(k)
1481 EXIT input_loop
1482 ELSE IF ( (z .GT. zprof(k)) .AND. &
1483 (z .LT. zprof(k+1)) ) THEN
1484 wgt0 = (z - zprof(k+1)) / &
1485 (zprof(k) - zprof(k+1))
1486 stor = MAX( wgt0 *cprof(k ) + &
1487 (1.-wgt0)*cprof(k+1), 0.)
1488 EXIT input_loop
1489 ENDIF
1490 ENDDO input_loop
1491 ENDIF
1492
1493 ! Here is where the chemistry value is constructed
1494 chem = fracref(nch-1)*stor*1.E6
1495
1496 ! special code for sulfate/h2so4
1497 if(nch.eq.p_sulf.and.p_nu0.gt.1)then
1498 chem=chem*(1.-so4vaptoaer)
1499 endif
1500
1501 RETURN
1502 END SUBROUTINE bdy_chem_value
1503 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1504
1505 #if (EM_CORE == 1 )
1506 SUBROUTINE flow_dep_bdy_chem ( chem, &
1507 chem_bxs,chem_btxs, &
1508 chem_bxe,chem_btxe, &
1509 chem_bys,chem_btys, &
1510 chem_bye,chem_btye, &
1511 dt, &
1512 spec_bdy_width,z, &
1513 have_bcs_chem, &
1514 u, v, config_flags, alt, &
1515 t,pb,p,t0,p1000mb,rcp,ph,phb,g, &
1516 spec_zone, ic, &
1517 ids,ide, jds,jde, kds,kde, & ! domain dims
1518 ims,ime, jms,jme, kms,kme, & ! memory dims
1519 ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1520 its,ite, jts,jte, kts,kte )
1521
1522 ! This subroutine sets zero gradient conditions for outflow and a set profile value
1523 ! for inflow in the boundary specified region. Note that field must be unstaggered.
1524 ! The velocities, u and v, will only be used to check their sign (coupled vels OK)
1525 ! spec_zone is the width of the outer specified b.c.s that are set here.
1526 ! (JD August 2000)
1527
1528 IMPLICIT NONE
1529
1530 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
1531 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
1532 INTEGER, INTENT(IN ) :: ips,ipe, jps,jpe, kps,kpe
1533 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
1534 INTEGER, INTENT(IN ) :: spec_zone,spec_bdy_width,ic
1535 REAL, INTENT(IN ) :: dt
1536
1537
1538 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: chem
1539 REAL, DIMENSION( jms:jme , kds:kde , spec_bdy_width), INTENT(IN ) :: chem_bxs, chem_bxe, chem_btxs, chem_btxe
1540 REAL, DIMENSION( ims:ime , kds:kde , spec_bdy_width), INTENT(IN ) :: chem_bys, chem_bye, chem_btys, chem_btye
1541 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: z
1542 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alt
1543 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: u
1544 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: v
1545 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
1546 INTENT(IN ) :: &
1547 ph,phb,t,pb,p
1548 real, INTENT (IN) :: g,rcp,t0,p1000mb
1549 TYPE( grid_config_rec_type ) config_flags
1550
1551 INTEGER :: i, j, k, numgas
1552 INTEGER :: ibs, ibe, jbs, jbe, itf, jtf, ktf
1553 INTEGER :: i_inner, j_inner
1554 INTEGER :: b_dist
1555 integer :: itestbc, i_bdy_method
1556 real tempfac,convfac
1557 real :: chem_bv_def
1558 logical :: have_bcs_chem
1559
1560 chem_bv_def = conmin
1561 numgas = get_last_gas(config_flags%chem_opt)
1562 itestbc=0
1563 if(p_nu0.gt.1)itestbc=1
1564 ibs = ids
1565 ibe = ide-1
1566 itf = min(ite,ide-1)
1567 jbs = jds
1568 jbe = jde-1
1569 jtf = min(jte,jde-1)
1570 ktf = kde-1
1571
1572 ! i_bdy_method determines which "bdy_chem_value" routine to use
1573 ! 1=radm2 or racm gas for p_so2 <= ic <= p_ho2
1574 ! 2=sorgam aerosol for p_so4aj <= ic <= p_corn
1575 ! 3=cbmz gas for p_hcl <= ic <= p_isopo2
1576 ! OR p_dms <= ic <= p_mtf
1577 ! 4=mosaic aerosol for p_so4_a01 <= ic <= p_num_a01
1578 ! OR p_so4_a02 <= ic <= p_num_a02
1579 ! OR ...
1580 ! 5=tracer mode
1581 ! 0=none for all other ic values
1582 ! (note: some cbmz packages use dms,...,mtf while others do not)
1583 ! (note: different mosaic packages use different number of sections)
1584 i_bdy_method = 0
1585 if ((ic .ge. p_so2) .and. (ic .le. p_ho2)) then
1586 i_bdy_method = 1
1587
1588 if (config_flags%chem_opt == RACM_KPP .or. &
1589 config_flags%chem_opt == RACMSORG_KPP .or. &
1590 config_flags%chem_opt == RACM_MIM_KPP .or. &
1591 config_flags%chem_opt == RACMSORG_KPP ) then
1592 i_bdy_method = 9
1593 end if
1594 if (config_flags%chem_opt == RACMPM_KPP ) then
1595 i_bdy_method = 9
1596 end if
1597
1598
1599 else if ((ic .ge. p_so4aj) .and. (ic .le. p_corn)) then
1600 i_bdy_method = 2
1601 else if ((ic .ge. p_hcl) .and. (ic .le. p_isopo2)) then
1602 i_bdy_method = 3
1603 else if ((ic .ge. p_dms) .and. (ic .le. p_mtf)) then
1604 i_bdy_method = 3
1605 else if ((ic .ge. p_so4_a01) .and. (ic .le. p_num_a01)) then
1606 i_bdy_method = 4
1607 else if ((ic .ge. p_so4_a02) .and. (ic .le. p_num_a02)) then
1608 i_bdy_method = 4
1609 else if ((ic .ge. p_so4_a03) .and. (ic .le. p_num_a03)) then
1610 i_bdy_method = 4
1611 else if ((ic .ge. p_so4_a04) .and. (ic .le. p_num_a04)) then
1612 i_bdy_method = 4
1613 else if ((ic .ge. p_so4_a05) .and. (ic .le. p_num_a05)) then
1614 i_bdy_method = 4
1615 else if ((ic .ge. p_so4_a06) .and. (ic .le. p_num_a06)) then
1616 i_bdy_method = 4
1617 else if ((ic .ge. p_so4_a07) .and. (ic .le. p_num_a07)) then
1618 i_bdy_method = 4
1619 else if ((ic .ge. p_so4_a08) .and. (ic .le. p_num_a08)) then
1620 i_bdy_method = 4
1621 else if (config_flags%chem_opt == CHEM_TRACER) then
1622 i_bdy_method = 5
1623 end if
1624 if (have_bcs_chem) i_bdy_method =6
1625 if (ic .lt. param_first_scalar) i_bdy_method = 0
1626
1627 !----------------------------------------------------------------------
1628 ! if (i_bdy_method .eq. 1) then
1629 ! print 90010, '_bdy_radm2 for ic=', ic, i_bdy_method
1630 ! else if (i_bdy_method .eq. 2) then
1631 ! print 90010, '_bdy_sorgam for ic=', ic, i_bdy_method
1632 ! else if (i_bdy_method .eq. 3) then
1633 ! print 90010, '_bdy_cbmz for ic=', ic, i_bdy_method
1634 ! else if (i_bdy_method .eq. 4) then
1635 ! print 90010, '_bdy_mosaic for ic=', ic, i_bdy_method
1636 ! else if (i_bdy_method .eq. 5) then
1637 ! print 90010, '_bdy_tracer for ic=', ic, i_bdy_method
1638 ! else
1639 ! print 90010, '_bdy_NONE** for ic=', ic, i_bdy_method
1640 ! end if
1641 !90010 format( a, 2(1x,i5) )
1642 !90020 format( a, 1p, 2e12.2 )
1643 !----------------------------------------------------------------------
1644
1645 IF (jts - jbs .lt. spec_zone) THEN
1646 ! Y-start boundary
1647 DO j = jts, min(jtf,jbs+spec_zone-1)
1648 b_dist = j - jbs
1649 DO k = kts, ktf
1650 DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist)
1651 i_inner = max(i,ibs+spec_zone)
1652 i_inner = min(i_inner,ibe-spec_zone)
1653 IF(v(i,k,j) .lt. 0.)THEN
1654 chem(i,k,j) = chem(i_inner,k,jbs+spec_zone)
1655 ELSE
1656 if (i_bdy_method .eq. 1) then
1657 CALL bdy_chem_value ( &
1658 chem(i,k,j), z(i,k,j), ic, numgas )
1659 else if (i_bdy_method .eq. 9) then
1660 CALL bdy_chem_value_racm( &
1661 chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
1662 else if (i_bdy_method .eq. 2) then
1663 tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
1664 convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
1665 CALL bdy_chem_value_sorgam ( &
1666 chem(i,k,j), z(i,k,j), ic, config_flags, &
1667 alt(i,k,j),convfac,g)
1668 else if (i_bdy_method .eq. 3) then
1669 CALL bdy_chem_value_cbmz ( &
1670 chem(i,k,j), z(i,k,j), ic, numgas )
1671 else if (i_bdy_method .eq. 4) then
1672 CALL bdy_chem_value_mosaic ( &
1673 chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
1674 else if (i_bdy_method .eq. 5) then
1675 CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
1676 else if (i_bdy_method .eq. 6) then
1677 CALL bdy_chem_value_gcm ( chem(i,k,j),chem_bys(i,k,1),chem_btys(i,k,1),dt,ic)
1678 else
1679 chem(i,k,j) = chem_bv_def
1680 endif
1681 ENDIF
1682 ENDDO
1683 ENDDO
1684 ENDDO
1685 ENDIF
1686 IF (jbe - jtf .lt. spec_zone) THEN
1687 ! Y-end boundary
1688 DO j = max(jts,jbe-spec_zone+1), jtf
1689 b_dist = jbe - j
1690 DO k = kts, ktf
1691 DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist)
1692 i_inner = max(i,ibs+spec_zone)
1693 i_inner = min(i_inner,ibe-spec_zone)
1694 IF(v(i,k,j+1) .gt. 0.)THEN
1695 chem(i,k,j) = chem(i_inner,k,jbe-spec_zone)
1696 ELSE
1697 if (i_bdy_method .eq. 1) then
1698 CALL bdy_chem_value ( &
1699 chem(i,k,j), z(i,k,j), ic, numgas )
1700 else if (i_bdy_method .eq. 9) then
1701 CALL bdy_chem_value_racm ( &
1702 chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
1703 else if (i_bdy_method .eq. 2) then
1704 tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
1705 convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
1706 CALL bdy_chem_value_sorgam ( &
1707 chem(i,k,j), z(i,k,j), ic, config_flags, &
1708 alt(i,k,j),convfac,g)
1709 else if (i_bdy_method .eq. 3) then
1710 CALL bdy_chem_value_cbmz ( &
1711 chem(i,k,j), z(i,k,j), ic, numgas )
1712 else if (i_bdy_method .eq. 4) then
1713 CALL bdy_chem_value_mosaic ( &
1714 chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
1715 else if (i_bdy_method .eq. 5) then
1716 CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
1717 else if (i_bdy_method .eq. 6) then
1718 CALL bdy_chem_value_gcm ( chem(i,k,j),chem_bye(i,k,1),chem_btye(i,k,1),dt,ic)
1719 else
1720 chem(i,k,j) = chem_bv_def
1721 endif
1722 ENDIF
1723 ENDDO
1724 ENDDO
1725 ENDDO
1726 ENDIF
1727
1728 IF (its - ibs .lt. spec_zone) THEN
1729 ! X-start boundary
1730 DO i = its, min(itf,ibs+spec_zone-1)
1731 b_dist = i - ibs
1732 DO k = kts, ktf
1733 DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1734 j_inner = max(j,jbs+spec_zone)
1735 j_inner = min(j_inner,jbe-spec_zone)
1736 IF(u(i,k,j) .lt. 0.)THEN
1737 chem(i,k,j) = chem(ibs+spec_zone,k,j_inner)
1738 ELSE
1739 if (i_bdy_method .eq. 1) then
1740 CALL bdy_chem_value ( &
1741 chem(i,k,j), z(i,k,j), ic, numgas )
1742 else if (i_bdy_method .eq. 9) then
1743 CALL bdy_chem_value_racm ( &
1744 chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
1745 else if (i_bdy_method .eq. 2) then
1746 tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
1747 convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
1748 CALL bdy_chem_value_sorgam ( &
1749 chem(i,k,j), z(i,k,j), ic, config_flags, &
1750 alt(i,k,j),convfac,g)
1751 else if (i_bdy_method .eq. 3) then
1752 CALL bdy_chem_value_cbmz ( &
1753 chem(i,k,j), z(i,k,j), ic, numgas )
1754 else if (i_bdy_method .eq. 4) then
1755 CALL bdy_chem_value_mosaic ( &
1756 chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
1757 else if (i_bdy_method .eq. 5) then
1758 CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
1759 else if (i_bdy_method .eq. 6) then
1760 CALL bdy_chem_value_gcm ( chem(i,k,j),chem_bxs(j,k,1),chem_btxs(j,k,1),dt,ic)
1761 else
1762 chem(i,k,j) = chem_bv_def
1763 endif
1764 ENDIF
1765 ENDDO
1766 ENDDO
1767 ENDDO
1768 ENDIF
1769
1770 IF (ibe - itf .lt. spec_zone) THEN
1771 ! X-end boundary
1772 DO i = max(its,ibe-spec_zone+1), itf
1773 b_dist = ibe - i
1774 DO k = kts, ktf
1775 DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
1776 j_inner = max(j,jbs+spec_zone)
1777 j_inner = min(j_inner,jbe-spec_zone)
1778 IF(u(i+1,k,j) .gt. 0.)THEN
1779 chem(i,k,j) = chem(ibe-spec_zone,k,j_inner)
1780 ELSE
1781 if (i_bdy_method .eq. 1) then
1782 CALL bdy_chem_value ( &
1783 chem(i,k,j), z(i,k,j), ic, numgas )
1784 else if (i_bdy_method .eq. 9) then
1785 CALL bdy_chem_value_racm ( &
1786 chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
1787 else if (i_bdy_method .eq. 2) then
1788 tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
1789 convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
1790 CALL bdy_chem_value_sorgam ( &
1791 chem(i,k,j), z(i,k,j), ic, config_flags, &
1792 alt(i,k,j),convfac,g)
1793 else if (i_bdy_method .eq. 3) then
1794 CALL bdy_chem_value_cbmz ( &
1795 chem(i,k,j), z(i,k,j), ic, numgas )
1796 else if (i_bdy_method .eq. 4) then
1797 CALL bdy_chem_value_mosaic ( &
1798 chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
1799 else if (i_bdy_method .eq. 5) then
1800 CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
1801 else if (i_bdy_method .eq. 6) then
1802 CALL bdy_chem_value_gcm ( chem(i,k,j),chem_bxe(j,k,1),chem_btxe(j,k,1),dt,ic)
1803 else
1804 chem(i,k,j) = chem_bv_def
1805 endif
1806 ENDIF
1807 ENDDO
1808 ENDDO
1809 ENDDO
1810 ENDIF
1811
1812 END SUBROUTINE flow_dep_bdy_chem
1813 #else
1814 SUBROUTINE flow_dep_bdy_chem ( chem, chem_b,chem_bt,dt, &
1815 spec_bdy_width,z, &
1816 ijds, ijde,have_bcs_chem, &
1817 u, v, config_flags, alt, &
1818 t,pb,p,t0,p1000mb,rcp,ph,phb,g, &
1819 spec_zone, ic, &
1820 ids,ide, jds,jde, kds,kde, & ! domain dims
1821 ims,ime, jms,jme, kms,kme, & ! memory dims
1822 ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1823 its,ite, jts,jte, kts,kte )
1824
1825 ! This subroutine sets zero gradient conditions for outflow and a set profile value
1826 ! for inflow in the boundary specified region. Note that field must be unstaggered.
1827 ! The velocities, u and v, will only be used to check their sign (coupled vels OK)
1828 ! spec_zone is the width of the outer specified b.c.s that are set here.
1829 ! (JD August 2000)
1830
1831 IMPLICIT NONE
1832
1833 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
1834 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
1835 INTEGER, INTENT(IN ) :: ips,ipe, jps,jpe, kps,kpe
1836 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
1837 INTEGER, INTENT(IN ) :: ijds,ijde
1838 INTEGER, INTENT(IN ) :: spec_zone,spec_bdy_width,ic
1839 REAL, INTENT(IN ) :: dt
1840
1841
1842 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: chem
1843 REAL, DIMENSION( ijds:ijde , kds:kde , spec_bdy_width, 4 ), INTENT(IN ) :: chem_b
1844 REAL, DIMENSION( ijds:ijde , kds:kde , spec_bdy_width, 4 ), INTENT(IN ) :: chem_bt
1845 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: z
1846 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alt
1847 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: u
1848 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: v
1849 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
1850 INTENT(IN ) :: &
1851 ph,phb,t,pb,p
1852 real, INTENT (IN) :: g,rcp,t0,p1000mb
1853 TYPE( grid_config_rec_type ) config_flags
1854
1855 INTEGER :: i, j, k, numgas
1856 INTEGER :: ibs, ibe, jbs, jbe, itf, jtf, ktf
1857 INTEGER :: i_inner, j_inner
1858 INTEGER :: b_dist
1859 integer :: itestbc, i_bdy_method
1860 real tempfac,convfac
1861 real :: chem_bv_def
1862 logical :: have_bcs_chem
1863
1864 chem_bv_def = conmin
1865 numgas = get_last_gas(config_flags%chem_opt)
1866 itestbc=0
1867 if(p_nu0.gt.1)itestbc=1
1868 ibs = ids
1869 ibe = ide-1
1870 itf = min(ite,ide-1)
1871 jbs = jds
1872 jbe = jde-1
1873 jtf = min(jte,jde-1)
1874 ktf = kde-1
1875
1876 ! i_bdy_method determines which "bdy_chem_value" routine to use
1877 ! 1=radm2 or racm gas for p_so2 <= ic <= p_ho2
1878 ! 2=sorgam aerosol for p_so4aj <= ic <= p_corn
1879 ! 3=cbmz gas for p_hcl <= ic <= p_isopo2
1880 ! OR p_dms <= ic <= p_mtf
1881 ! 4=mosaic aerosol for p_so4_a01 <= ic <= p_num_a01
1882 ! OR p_so4_a02 <= ic <= p_num_a02
1883 ! OR ...
1884 ! 5=tracer mode
1885 ! 0=none for all other ic values
1886 ! (note: some cbmz packages use dms,...,mtf while others do not)
1887 ! (note: different mosaic packages use different number of sections)
1888 i_bdy_method = 0
1889 if ((ic .ge. p_so2) .and. (ic .le. p_ho2)) then
1890 i_bdy_method = 1
1891
1892 if (config_flags%chem_opt == RACM_KPP .or. &
1893 config_flags%chem_opt == RACMSORG_KPP .or. &
1894 config_flags%chem_opt == RACM_MIM_KPP .or. &
1895 config_flags%chem_opt == RACMSORG_KPP ) then
1896 i_bdy_method = 9
1897 end if
1898 if (config_flags%chem_opt == RACMPM_KPP ) then
1899 i_bdy_method = 9
1900 end if
1901
1902
1903 else if ((ic .ge. p_so4aj) .and. (ic .le. p_corn)) then
1904 i_bdy_method = 2
1905 else if ((ic .ge. p_hcl) .and. (ic .le. p_isopo2)) then
1906 i_bdy_method = 3
1907 else if ((ic .ge. p_dms) .and. (ic .le. p_mtf)) then
1908 i_bdy_method = 3
1909 else if ((ic .ge. p_so4_a01) .and. (ic .le. p_num_a01)) then
1910 i_bdy_method = 4
1911 else if ((ic .ge. p_so4_a02) .and. (ic .le. p_num_a02)) then
1912 i_bdy_method = 4
1913 else if ((ic .ge. p_so4_a03) .and. (ic .le. p_num_a03)) then
1914 i_bdy_method = 4
1915 else if ((ic .ge. p_so4_a04) .and. (ic .le. p_num_a04)) then
1916 i_bdy_method = 4
1917 else if ((ic .ge. p_so4_a05) .and. (ic .le. p_num_a05)) then
1918 i_bdy_method = 4
1919 else if ((ic .ge. p_so4_a06) .and. (ic .le. p_num_a06)) then
1920 i_bdy_method = 4
1921 else if ((ic .ge. p_so4_a07) .and. (ic .le. p_num_a07)) then
1922 i_bdy_method = 4
1923 else if ((ic .ge. p_so4_a08) .and. (ic .le. p_num_a08)) then
1924 i_bdy_method = 4
1925 else if (config_flags%chem_opt == CHEM_TRACER) then
1926 i_bdy_method = 5
1927 end if
1928 if (have_bcs_chem) i_bdy_method =6
1929 if (ic .lt. param_first_scalar) i_bdy_method = 0
1930
1931 !----------------------------------------------------------------------
1932 ! if (i_bdy_method .eq. 1) then
1933 ! print 90010, '_bdy_radm2 for ic=', ic, i_bdy_method
1934 ! else if (i_bdy_method .eq. 2) then
1935 ! print 90010, '_bdy_sorgam for ic=', ic, i_bdy_method
1936 ! else if (i_bdy_method .eq. 3) then
1937 ! print 90010, '_bdy_cbmz for ic=', ic, i_bdy_method
1938 ! else if (i_bdy_method .eq. 4) then
1939 ! print 90010, '_bdy_mosaic for ic=', ic, i_bdy_method
1940 ! else if (i_bdy_method .eq. 5) then
1941 ! print 90010, '_bdy_tracer for ic=', ic, i_bdy_method
1942 ! else
1943 ! print 90010, '_bdy_NONE** for ic=', ic, i_bdy_method
1944 ! end if
1945 !90010 format( a, 2(1x,i5) )
1946 !90020 format( a, 1p, 2e12.2 )
1947 !----------------------------------------------------------------------
1948
1949 ! if(ic.eq.p_O3)THEN
1950 ! write(0,*)'in flow_chem ',jts,jbs,spec_zone
1951 ! write(0,*)'in flow_chem ',its,ibs,b_dist,i_bdy_method,ic
1952 ! endif
1953 IF (jts - jbs .lt. spec_zone) THEN
1954 ! Y-start boundary
1955 DO j = jts, min(jtf,jbs+spec_zone-1)
1956 b_dist = j - jbs
1957 DO k = kts, ktf
1958 DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist)
1959 i_inner = max(i,ibs+spec_zone)
1960 i_inner = min(i_inner,ibe-spec_zone)
1961 IF(v(i,k,j) .lt. 0.)THEN
1962 chem(i,k,j) = chem(i_inner,k,jbs+spec_zone)
1963 ! if(j.eq.jts+1.and.k.eq.kts.and.ic.eq.p_o3)then
1964 ! write(0,*)'Yflow',i,j,k,i_bdy_method
1965 ! write(0,*)chem(i_inner,k,jbs+spec_zone),v(i,k,j)
1966 ! endif
1967 ELSE
1968 if (i_bdy_method .eq. 1) then
1969 CALL bdy_chem_value ( &
1970 chem(i,k,j), z(i,k,j), ic, numgas )
1971 else if (i_bdy_method .eq. 9) then
1972 CALL bdy_chem_value_racm( &
1973 chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
1974 else if (i_bdy_method .eq. 2) then
1975 tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
1976 convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
1977 CALL bdy_chem_value_sorgam ( &
1978 chem(i,k,j), z(i,k,j), ic, config_flags, &
1979 alt(i,k,j),convfac,g)
1980 else if (i_bdy_method .eq. 3) then
1981 CALL bdy_chem_value_cbmz ( &
1982 chem(i,k,j), z(i,k,j), ic, numgas )
1983 else if (i_bdy_method .eq. 4) then
1984 CALL bdy_chem_value_mosaic ( &
1985 chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
1986 else if (i_bdy_method .eq. 5) then
1987 CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
1988 else if (i_bdy_method .eq. 6) then
1989 CALL bdy_chem_value_gcm ( chem(i,k,j),chem_b(i,k,1,P_YSB),chem_bt(i,k,1,P_YSB),dt,ic)
1990 ! if(k.eq.kts.and.ic.eq.p_o3)then
1991 ! write(0,*)'Ygcm',i,j,k,i_bdy_method
1992 ! write(0,*)chem(i,k,j),chem_b(i,k,1,P_YSB),chem_bt(i,k,1,P_YSB),dt
1993 ! endif
1994 else
1995 chem(i,k,j) = chem_bv_def
1996 endif
1997 ENDIF
1998 ENDDO
1999 ENDDO
2000 ENDDO
2001 ENDIF
2002 IF (jbe - jtf .lt. spec_zone) THEN
2003 ! Y-end boundary
2004 DO j = max(jts,jbe-spec_zone+1), jtf
2005 b_dist = jbe - j
2006 DO k = kts, ktf
2007 DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist)
2008 i_inner = max(i,ibs+spec_zone)
2009 i_inner = min(i_inner,ibe-spec_zone)
2010 IF(v(i,k,j+1) .gt. 0.)THEN
2011 chem(i,k,j) = chem(i_inner,k,jbe-spec_zone)
2012 ELSE
2013 if (i_bdy_method .eq. 1) then
2014 CALL bdy_chem_value ( &
2015 chem(i,k,j), z(i,k,j), ic, numgas )
2016 else if (i_bdy_method .eq. 9) then
2017 CALL bdy_chem_value_racm ( &
2018 chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2019 else if (i_bdy_method .eq. 2) then
2020 tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2021 convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2022 CALL bdy_chem_value_sorgam ( &
2023 chem(i,k,j), z(i,k,j), ic, config_flags, &
2024 alt(i,k,j),convfac,g)
2025 else if (i_bdy_method .eq. 3) then
2026 CALL bdy_chem_value_cbmz ( &
2027 chem(i,k,j), z(i,k,j), ic, numgas )
2028 else if (i_bdy_method .eq. 4) then
2029 CALL bdy_chem_value_mosaic ( &
2030 chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2031 else if (i_bdy_method .eq. 5) then
2032 CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2033 else if (i_bdy_method .eq. 6) then
2034 CALL bdy_chem_value_gcm ( chem(i,k,j),chem_b(i,k,1,P_YEB),chem_bt(i,k,1,P_YEB),dt,ic)
2035 else
2036 chem(i,k,j) = chem_bv_def
2037 endif
2038 ENDIF
2039 ENDDO
2040 ENDDO
2041 ENDDO
2042 ENDIF
2043
2044 IF (its - ibs .lt. spec_zone) THEN
2045 ! X-start boundary
2046 DO i = its, min(itf,ibs+spec_zone-1)
2047 b_dist = i - ibs
2048 DO k = kts, ktf
2049 DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2050 j_inner = max(j,jbs+spec_zone)
2051 j_inner = min(j_inner,jbe-spec_zone)
2052 IF(u(i,k,j) .lt. 0.)THEN
2053 chem(i,k,j) = chem(ibs+spec_zone,k,j_inner)
2054 ELSE
2055 if (i_bdy_method .eq. 1) then
2056 CALL bdy_chem_value ( &
2057 chem(i,k,j), z(i,k,j), ic, numgas )
2058 else if (i_bdy_method .eq. 9) then
2059 CALL bdy_chem_value_racm ( &
2060 chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2061 else if (i_bdy_method .eq. 2) then
2062 tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2063 convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2064 CALL bdy_chem_value_sorgam ( &
2065 chem(i,k,j), z(i,k,j), ic, config_flags, &
2066 alt(i,k,j),convfac,g)
2067 else if (i_bdy_method .eq. 3) then
2068 CALL bdy_chem_value_cbmz ( &
2069 chem(i,k,j), z(i,k,j), ic, numgas )
2070 else if (i_bdy_method .eq. 4) then
2071 CALL bdy_chem_value_mosaic ( &
2072 chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2073 else if (i_bdy_method .eq. 5) then
2074 CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2075 else if (i_bdy_method .eq. 6) then
2076 CALL bdy_chem_value_gcm ( chem(i,k,j),chem_b(j,k,1,P_XSB),chem_bt(j,k,1,P_XSB),dt,ic)
2077 else
2078 chem(i,k,j) = chem_bv_def
2079 endif
2080 ENDIF
2081 ENDDO
2082 ENDDO
2083 ENDDO
2084 ENDIF
2085
2086 IF (ibe - itf .lt. spec_zone) THEN
2087 ! X-end boundary
2088 DO i = max(its,ibe-spec_zone+1), itf
2089 b_dist = ibe - i
2090 DO k = kts, ktf
2091 DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2092 j_inner = max(j,jbs+spec_zone)
2093 j_inner = min(j_inner,jbe-spec_zone)
2094 IF(u(i+1,k,j) .gt. 0.)THEN
2095 chem(i,k,j) = chem(ibe-spec_zone,k,j_inner)
2096 ELSE
2097 if (i_bdy_method .eq. 1) then
2098 CALL bdy_chem_value ( &
2099 chem(i,k,j), z(i,k,j), ic, numgas )
2100 else if (i_bdy_method .eq. 9) then
2101 CALL bdy_chem_value_racm ( &
2102 chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2103 else if (i_bdy_method .eq. 2) then
2104 tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2105 convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2106 CALL bdy_chem_value_sorgam ( &
2107 chem(i,k,j), z(i,k,j), ic, config_flags, &
2108 alt(i,k,j),convfac,g)
2109 else if (i_bdy_method .eq. 3) then
2110 CALL bdy_chem_value_cbmz ( &
2111 chem(i,k,j), z(i,k,j), ic, numgas )
2112 else if (i_bdy_method .eq. 4) then
2113 CALL bdy_chem_value_mosaic ( &
2114 chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2115 else if (i_bdy_method .eq. 5) then
2116 CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2117 else if (i_bdy_method .eq. 6) then
2118 CALL bdy_chem_value_gcm ( chem(i,k,j),chem_b(j,k,1,P_XEB),chem_bt(j,k,1,P_XEB),dt,ic)
2119 else
2120 chem(i,k,j) = chem_bv_def
2121 endif
2122 ENDIF
2123 ENDDO
2124 ENDDO
2125 ENDDO
2126 ENDIF
2127
2128 END SUBROUTINE flow_dep_bdy_chem
2129 #endif
2130 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2131 SUBROUTINE bdy_chem_value_gcm ( chem, chem_b, chem_bt, dt,ic)
2132
2133 IMPLICIT NONE
2134
2135 REAL, intent(OUT) :: chem
2136 REAL, intent(IN) :: chem_b
2137 REAL, intent(IN) :: chem_bt
2138 REAL, intent(IN) :: dt
2139 INTEGER, intent(IN) :: ic
2140
2141
2142 CHARACTER (LEN=80) :: message
2143 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2144
2145 ! if( nch .GT. numgas) then
2146 ! message = ' Input_chem_profile: wrong number of chemical species'
2147 ! return
2148 ! CALL WRF_ERROR_FATAL ( message )
2149 ! endif
2150
2151
2152 !print*,'before',chem,chem_bt ,dt, ic
2153
2154 chem=max(epsilc,chem_b + chem_bt * dt)
2155 !print*,'after',chem
2156 RETURN
2157 END SUBROUTINE bdy_chem_value_gcm
2158 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2159
2160 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2161 SUBROUTINE cv_mmdd_jday ( YY, MM, DD, JDAY)
2162 !
2163 ! Subroutine to compute the julian day given the month and day
2164 !
2165 !
2166 INTEGER, INTENT(IN ) :: YY, MM, DD
2167 INTEGER, INTENT(OUT) :: JDAY
2168
2169 INTEGER, DIMENSION(12) :: imon, imon_a
2170 INTEGER :: i
2171
2172 DATA imon_a /0,31,59,90,120,151,181,212,243,273,304,334/
2173 !
2174 !..... Check for leap year.
2175 !
2176 do i=1,12
2177 imon(i) = imon_a(i)
2178 enddo
2179 if(YY .eq. (YY/4)*4) then
2180 do i=3,12
2181 imon(i) = imon(i) + 1
2182 enddo
2183 endif
2184 !
2185 !..... Convert month, day to julian day.
2186 !
2187 jday = imon(mm) + dd
2188
2189
2190 END SUBROUTINE cv_mmdd_jday
2191
2192 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2193
2194 integer FUNCTION get_last_gas(chem_opt)
2195 implicit none
2196 integer, intent(in) :: chem_opt
2197
2198 ! Determine the index of the last gas species, which depends
2199 ! upon the gas mechanism.
2200
2201 select case (chem_opt)
2202 case (0)
2203 get_last_gas = 0
2204 case (RADM2,RADM2_KPP,RADM2SORG,RADM2SORG_KPP,RACM,RACMSORG,RACM_KPP,RACMPM_KPP,RACMSORG_KPP, RACM_MIM_KPP, RADM2SORG_AQ,RACMSORG_AQ)
2205 get_last_gas = p_ho2
2206
2207 case (CBMZ)
2208 get_last_gas = p_mtf
2209
2210 case (CBMZ_BB,CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN,CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ)
2211 get_last_gas = p_isopo2
2212
2213 case (CHEM_TRACER)
2214 get_last_gas = p_co
2215
2216 case default
2217 call wrf_error_fatal("get_last_gas: could not decipher chem_opt value")
2218
2219 end select
2220
2221 END FUNCTION get_last_gas
2222 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2223
2224 !****************************************************************
2225 ! *
2226 ! SUBROUTINE TO SET AEROSOL BC VALUES USING THE *
2227 ! aer_bc_opt == aer_bc_pnnl OPTION. *
2228 ! *
2229 ! wig 22-Apr-2004, original routine *
2230 ! rce 25-apr-2004 - changed name to *
2231 ! "sorgam_set_aer_bc_pnnl" *
2232 ! wig 7-May-2004, added height dependance *
2233 ! *
2234 ! CALLS THE FOLLOWING SUBROUTINES: NONE *
2235 ! *
2236 ! CALLED BY: bdy_chem_value_sorgam *
2237 ! *
2238 !****************************************************************
2239 SUBROUTINE sorgam_set_aer_bc_pnnl( chem, z, nch )
2240 USE module_data_sorgam
2241 implicit none
2242
2243 INTEGER,INTENT(IN ) :: nch
2244 real,intent(in ) :: z
2245 REAL,INTENT(INOUT ) :: chem
2246
2247 REAL :: mult, &
2248 m3acc, m3cor, m3nuc, &
2249 bv_so4ai, bv_so4aj, &
2250 bv_nh4ai, bv_nh4aj, &
2251 bv_no3ai, bv_no3aj, &
2252 bv_eci, bv_ecj, &
2253 bv_p25i, bv_p25j, &
2254 bv_orgpai,bv_orgpaj, &
2255 bv_antha, bv_seas, bv_soila
2256
2257 !
2258 ! Determine height multiplier...
2259 ! This should mimic the calculation in sorgam_init_aer_ic_pnnl,
2260 ! mosaic_init_wrf_mixrats_opt2, and bdy_chem_value_mosaic
2261 !!$! Height(m) Multiplier
2262 !!$! --------- ----------
2263 !!$! <=2000 1.0
2264 !!$! 2000<z<3000 linear transition zone to 0.5
2265 !!$! 3000<z<5000 linear transision zone to 0.25
2266 !!$! >=5000 0.25
2267 !!$!
2268 !!$! which translates to:
2269 !!$! 2000<z<3000 mult = 1.0 + (z-2000.)*(0.5-1.0)/(3000.-2000.)
2270 !!$! 3000<z<5000 mult = 0.5 + (z-3000.)*(0.25-0.5)/(5000.-3000.)
2271 !!$! or in reduced form:
2272 !!$ if( z <= 2000. ) then
2273 !!$ mult = 1.0
2274 !!$ elseif( z > 2000. &
2275 !!$ .and. z <= 3000. ) then
2276 !!$ mult = 1.0 - 0.0005*(z-2000.)
2277 !!$ elseif( z > 3000. &
2278 !!$ .and. z <= 5000. ) then
2279 !!$ mult = 0.5 - 1.25e-4*(z-3000.)
2280 !!$ else
2281 !!$ mult = 0.25
2282 !!$ end if
2283 ! Updated aerosol profile multiplier 1-Apr-2005:
2284 ! Height(m) Multiplier
2285 ! --------- ----------
2286 ! <=2000 1.0
2287 ! 2000<z<3000 linear transition zone to 0.25
2288 ! 3000<z<5000 linear transision zone to 0.125
2289 ! >=5000 0.125
2290 !
2291 ! which translates to:
2292 ! 2000<z<3000 mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
2293 ! 3000<z<5000 mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
2294 ! or in reduced form:
2295 if( z <= 2000. ) then
2296 mult = 1.0
2297 elseif( z > 2000. &
2298 .and. z <= 3000. ) then
2299 mult = 1.0 - 0.00075*(z-2000.)
2300 elseif( z > 3000. &
2301 .and. z <= 5000. ) then
2302 mult = 0.25 - 4.166666667e-5*(z-3000.)
2303 else
2304 mult = 0.125
2305 end if
2306
2307 ! These should match what is in sorgam_init_aer_ic_pnnl.
2308 ! Boundary values as of 2-Dec-2004:
2309 bv_so4aj = mult*2.375
2310 bv_so4ai = mult*0.179
2311 bv_nh4aj = mult*0.9604
2312 bv_nh4ai = mult*0.0196
2313 bv_no3aj = mult*0.0650
2314 bv_no3ai = mult*0.0050
2315 bv_ecj = mult*0.1630
2316 bv_eci = mult*0.0120
2317 bv_p25j = mult*0.6350
2318 bv_p25i = mult*0.0490
2319 bv_orgpaj = mult*0.9300
2320 bv_orgpai = mult*0.0700
2321 bv_antha = mult*2.2970
2322 bv_seas = mult*0.2290
2323 bv_soila = conmin
2324
2325 ! m3... calculations should match the very end of module_aerosols_sorgam.F
2326 !... i-mode (note that the 8 SOA species have bv=conmin)
2327 m3nuc = so4fac*bv_so4ai + nh4fac*bv_nh4ai + &
2328 no3fac*bv_no3ai + &
2329 orgfac*8.0*conmin + orgfac*bv_orgpai + &
2330 anthfac*bv_p25i + anthfac*bv_eci
2331
2332 !... j-mode (note that the 8 SOA species have bv=conmin)
2333 m3acc = so4fac*bv_so4aj + nh4fac*bv_nh4aj + &
2334 no3fac*bv_no3aj + &
2335 orgfac*8.0*conmin + orgfac*bv_orgpaj + &
2336 anthfac*bv_p25j + anthfac*bv_ecj
2337
2338 !...c-mode
2339 m3cor = soilfac*bv_soila + seasfac*bv_seas + &
2340 anthfac*bv_antha
2341
2342 ! Cannot set_sulf here because it is a "radm2" species whose bc value
2343 ! is set via bdy_chem_value. Instead, xl(iref(p_sulf-1),:) is set to
2344 ! the value conmin in subroutine gasprofile_init_pnnl
2345 ! if( nch == p_sulf ) chem = conmin !as per rce's 0 recommendation
2346
2347 if( nch == p_so4aj ) chem = bv_so4aj
2348 if( nch == p_so4ai ) chem = bv_so4ai
2349 if( nch == p_nh4aj ) chem = bv_nh4aj
2350 if( nch == p_nh4ai ) chem = bv_nh4ai
2351 if( nch == p_no3aj ) chem = bv_no3aj
2352 if( nch == p_no3ai ) chem = bv_no3ai
2353 if( nch == p_ecj ) chem = bv_ecj
2354 if( nch == p_eci ) chem = bv_eci
2355 if( nch == p_p25j ) chem = bv_p25j
2356 if( nch == p_p25i ) chem = bv_p25i
2357 if( nch == p_orgpaj ) chem = bv_orgpaj
2358 if( nch == p_orgpai ) chem = bv_orgpai
2359
2360 if( nch == p_orgaro1j) chem = conmin
2361 if( nch == p_orgaro1i) chem = conmin
2362 if( nch == p_orgaro2j) chem = conmin
2363 if( nch == p_orgaro2i) chem = conmin
2364 if( nch == p_orgalk1j) chem = conmin
2365 if( nch == p_orgalk1i) chem = conmin
2366 if( nch == p_orgole1j) chem = conmin
2367 if( nch == p_orgole1i) chem = conmin
2368 if( nch == p_orgba1j ) chem = conmin
2369 if( nch == p_orgba1i ) chem = conmin
2370 if( nch == p_orgba2j ) chem = conmin
2371 if( nch == p_orgba2i ) chem = conmin
2372 if( nch == p_orgba3j ) chem = conmin
2373 if( nch == p_orgba3i ) chem = conmin
2374 if( nch == p_orgba4j ) chem = conmin
2375 if( nch == p_orgba4i ) chem = conmin
2376
2377 if( nch == p_antha ) chem = bv_antha
2378 if( nch == p_soila ) chem = bv_soila
2379 if( nch == p_seas ) chem = bv_seas
2380
2381 if( nch == p_nu0 ) chem = m3nuc/((dginin**3)*esn36)
2382 if( nch == p_ac0 ) chem = m3acc/((dginia**3)*esa36)
2383 if( nch == p_corn ) chem = m3cor/((dginic**3)*esc36)
2384
2385 END SUBROUTINE sorgam_set_aer_bc_pnnl
2386 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2387
2388 !****************************************************************
2389 ! *
2390 ! SUBROUTINE TO OVERWRITE THE PREDEFINED OZONE PROFILE *
2391 ! WHEN gas_ic_opt == gas_ic_pnnl *
2392 ! OR gas_bc_opt == gas_bc_pnnl *
2393 ! *
2394 ! wig, 21-Apr-2004 *
2395 ! rce 25-apr-2004 - changed name to *
2396 ! "gasprofile_init_pnnl" *
2397 ! *
2398 ! CALLS THE FOLLOWING SUBROUTINES: NONE *
2399 ! *
2400 ! CALLED BY: chem_init *
2401 ! input_chem_profile *
2402 ! *
2403 !****************************************************************
2404 SUBROUTINE gasprofile_init_pnnl
2405 use module_data_sorgam,only: conmin
2406 implicit none
2407 integer :: k
2408
2409 call wrf_debug ( 500 , 'wrfchem:gasprofile_init_pnnl' )
2410 ! print*,'gasprofile_init_pnnl redefining o3 and sulf profiles.'
2411
2412 ! Original O3 profile values:
2413 ! / 1.68E-07, 1.68E-07, 5.79E-08, 5.24E-08, 5.26E-08, &
2414 ! 5.16E-08, 4.83E-08, 4.50E-08, 4.16E-08, 3.80E-08, 3.56E-08, &
2415 ! 3.35E-08, 3.15E-08, 3.08E-08, 3.06E-08, 3.00E-08/
2416
2417 ! Note that heights associated with 2nd index of xl correspond to upside-down
2418 ! zfa values that have been "de-staggered".
2419 ! Height = 0.5*(zfa(1:kx) + zfa(2:kx+1)) and then flipped:
2420 ! / 18500., 14050., 11150., 9355., 7705., 6120., 4675., 3430.,
2421 ! 2430., 1720., 1195., 781.5, 494., 298.5, 148.5, 42.5
2422
2423 if( p_o3 > 1 ) then
2424 #if (CASENAME == 0)
2425 !Rounded to closest level:
2426 xl(iref(p_o3-1),11:16) = 4.00e-8 !40 ppbv below 1 km
2427 xl(iref(p_o3-1),3:10) = 6.50e-8 !65 ppbv > 2 km and < stratosphere
2428 ! Changed from 70 ppbv 1-Apr-2005
2429 #endif
2430 #if (CASENAME == 1)
2431 xl(iref(p_o3-1),11:16) = 3.50e-8 !35 ppbv below 1 km
2432 xl(iref(p_o3-1),3:10) = 6.00e-8 !60 ppbv > 2 km and < stratosphere
2433 #endif
2434 end if
2435
2436 #if (CASENAME == 1)
2437 ! so2 profile based on mirage 2 output, used for neaqs case, 7-20-05 egc
2438 ! decreased by one magnitude, 27-oct-2005 wig
2439 if( p_so2 > 1 ) then
2440 xl(iref(p_so2-1), 1:2) = 0.035e-10
2441 xl(iref(p_so2-1), 3) = 0.081e-10
2442 xl(iref(p_so2-1), 4:8) = 0.10e-10
2443 xl(iref(p_so2-1), 9) = 0.60e-10
2444 xl(iref(p_so2-1), 10) = 1.1e-10
2445 xl(iref(p_so2-1), 11) = 1.46e-10
2446 xl(iref(p_so2-1), 12) = 1.74e-10
2447 xl(iref(p_so2-1), 13) = 1.94e-10
2448 xl(iref(p_so2-1), 14) = 2.80e-10
2449 xl(iref(p_so2-1), 15:16) = 3.0e-10
2450 end if
2451 #endif
2452
2453 if( p_sulf > 1 ) then
2454 xl(iref(p_sulf-1),:) = conmin
2455 end if
2456
2457 end SUBROUTINE gasprofile_init_pnnl
2458
2459 #ifdef CHEM_DBG_I
2460 !-----------------------------------------------------------------------
2461 subroutine chem_dbg(i,j,k,dtstep,itimestep, &
2462 dz8w,t_phy,p_phy,rho_phy,chem, &
2463 e_so2,e_no,e_co,e_eth,e_hc3,e_hc5,e_hc8,e_xyl,e_ol2,e_olt, &
2464 e_oli,e_tol,e_csl,e_hcho,e_ald,e_ket,e_ora2,e_nh3, &
2465 e_pm10,e_pm25,e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
2466 e_no2,e_ch3oh,e_c2h5oh,e_iso, &
2467 e_so4j,e_so4c,e_no3j,e_no3c,e_orgc,e_ecc, &
2468 ids,ide, jds,jde, kds,kde, &
2469 ims,ime, jms,jme, kms,kme, &
2470 its,ite, jts,jte, kts,kte, &
2471 kemit, &
2472 ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
2473 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
2474 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
2475 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,ph_n2o5, &
2476 ph_o2 )
2477
2478 IMPLICIT NONE
2479 INTEGER, INTENT(IN ) :: i,j,k, &
2480 ids,ide, jds,jde, kds,kde, &
2481 ims,ime, jms,jme, kms,kme, &
2482 its,ite, jts,jte, kts,kte, &
2483 kemit
2484 real, intent(in ) :: dtstep
2485 integer, intent(in ) :: itimestep
2486 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
2487 INTENT(INOUT ) :: chem
2488 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
2489 INTENT(IN ) :: dz8w,t_phy,p_phy,rho_phy
2490 REAL, DIMENSION( ims:ime, kms:kemit, jms:jme ), &
2491 INTENT(IN ) :: &
2492 e_pm25i,e_pm25j,e_eci,e_ecj,e_orgi,e_orgj, &
2493 e_so2,e_no,e_co,e_eth,e_hc3,e_hc5,e_hc8,e_xyl,e_ol2,e_olt, &
2494 e_oli,e_tol,e_csl,e_hcho,e_ald,e_ket,e_ora2,e_pm25,e_pm10,e_nh3, &
2495 e_no2,e_ch3oh,e_c2h5oh,e_iso, &
2496 e_so4j,e_so4c,e_no3j,e_no3c,e_orgc,e_ecc
2497 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
2498 INTENT(IN ), OPTIONAL :: &
2499 ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
2500 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
2501 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
2502 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,ph_n2o5, &
2503 ph_o2
2504
2505 integer :: n
2506 real :: conva,convg
2507
2508 print*,"itimestep =",itimestep
2509
2510 print*,"MET DATA AT (i,k,j):",i,k,j
2511 print*,"t_phy,p_phy,rho_phy=",t_phy(i,k,j),p_phy(i,k,j),rho_phy(i,k,j)
2512
2513 if(dz8w(i,k,j) /= 0.) then
2514 conva = dtstep/(dz8w(i,k,j)*60.)
2515 convg = 4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.)
2516 print*,"ADJUSTED EMISSIONS (PPM) AT (i,k,j):",i,k,j
2517 print*,"dtstep,dz8w(i,k,j):",dtstep,dz8w(i,k,j)
2518 print*,"e_pm25 i,j:",e_pm25i(i,k,j)*conva, &
2519 e_pm25j(i,k,j)*conva
2520 print*,"e_ec i,j:",e_eci(i,k,j)*conva, &
2521 e_ecj(i,k,j)*conva
2522 print*,"e_org i,j:",e_orgi(i,k,j)*conva, &
2523 e_orgj(i,k,j)*conva
2524 print*,"e_so2:",e_so2(i,k,j)*convg
2525 print*,"e_no:",e_no(i,k,j)*convg
2526 print*,"e_co:",e_co(i,k,j)*convg
2527 print*,"e_eth:",e_eth(i,k,j)*convg
2528 print*,"e_hc3:",e_hc3(i,k,j)*convg
2529 print*,"e_hc5:",e_hc5(i,k,j)*convg
2530 print*,"e_hc8:",e_hc8(i,k,j)*convg
2531 print*,"e_xyl:",e_xyl(i,k,j)*convg
2532 print*,"e_ol2:",e_ol2(i,k,j)*convg
2533 print*,"e_olt:",e_olt(i,k,j)*convg
2534 print*,"e_oli:",e_oli(i,k,j)*convg
2535 print*,"e_tol:",e_tol(i,k,j)*convg
2536 print*,"e_csl:",e_csl(i,k,j)*convg
2537 print*,"e_hcho:",e_hcho(i,k,j)*convg
2538 print*,"e_ald:",e_ald(i,k,j)*convg
2539 print*,"e_ket:",e_ket(i,k,j)*convg
2540 print*,"e_ora2:",e_ora2(i,k,j)*convg
2541 print*,"e_pm25:",e_pm25(i,k,j)*conva
2542 print*,"e_pm10:",e_pm10(i,k,j)*conva
2543 print*,"e_nh3:",e_nh3(i,k,j)*convg
2544 print*,"e_no2:",e_no2(i,k,j)*convg
2545 print*,"e_ch3oh:",e_ch3oh(i,k,j)*convg
2546 print*,"e_c2h5oh:",e_c2h5oh(i,k,j)*convg
2547 print*,"e_iso:",e_iso(i,k,j)*convg
2548 print*,"e_so4 f,c:",e_so4j(i,k,j)*conva, &
2549 e_so4c(i,k,j)*conva
2550 print*,"e_no3 f,c:",e_no3j(i,k,j)*conva, &
2551 e_no3c(i,k,j)*conva
2552 print*,"e_orgc:",e_orgc(i,k,j)*conva
2553 print*,"e_ecc:",e_ecc(i,k,j)*conva
2554 print*
2555 else
2556 print*,"dz8w=0 so cannot show adjusted emissions"
2557 end if
2558 print*,"CHEM_DBG PRINT (PPM or ug/m^3) AT (i,k,j):",i,k,j
2559 do n=1,num_chem
2560 print*,n,chem(i,k,j,n)
2561 end do
2562 if( present(ph_macr) ) then
2563 print*,"PHOTOLYSIS DATA:"
2564 print*,"ph_macr:",ph_macr(i,:,j)
2565 print*,"ph_o31d:",ph_o31d(i,:,j)
2566 print*,"ph_o33p:",ph_o33p(i,:,j)
2567 print*,"ph_no2:",ph_no2(i,:,j)
2568 print*,"ph_no3o2:",ph_no3o2(i,:,j)
2569 print*,"ph_no3o:",ph_no3o(i,:,j)
2570 print*,"ph_hno2:",ph_hno2(i,:,j)
2571 print*,"ph_hno3:",ph_hno3(i,:,j)
2572 print*,"ph_hno4:",ph_hno4(i,:,j)
2573 print*,"ph_h2o2:",ph_h2o2(i,:,j)
2574 print*,"ph_ch2or:",ph_ch2or(i,:,j)
2575 print*,"ph_ch2om:",ph_ch2om(i,:,j)
2576 print*,"ph_ch3cho:",ph_ch3cho(i,:,j)
2577 print*,"ph_ch3coch3:",ph_ch3coch3(i,:,j)
2578 print*,"ph_ch3coc2h5:",ph_ch3coc2h5(i,:,j)
2579 print*,"ph_hcocho:",ph_hcocho(i,:,j)
2580 print*,"ph_ch3cocho:",ph_ch3cocho(i,:,j)
2581 print*,"ph_hcochest:",ph_hcochest(i,:,j)
2582 print*,"ph_ch3o2h:",ph_ch3o2h(i,:,j)
2583 print*,"ph_ch3coo2h:",ph_ch3coo2h(i,:,j)
2584 print*,"ph_ch3ono2:",ph_ch3ono2(i,:,j)
2585 print*,"ph_hcochob:",ph_hcochob(i,:,j)
2586 print*,"ph_n2o5:",ph_n2o5(i,:,j)
2587 print*,"ph_o2:",ph_o2(i,:,j)
2588 end if
2589 print*
2590 end subroutine chem_dbg
2591 #endif
2592
2593 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2594 SUBROUTINE med_read_bin_chem_emiss ( grid , config_flags ,intime, itime_max)
2595 ! Driver layer
2596 ! USE module_domain
2597 ! USE module_io_domain
2598 ! USE module_timing
2599 ! Model layer
2600 ! USE module_configure
2601 USE module_bc_time_utilities
2602 #ifdef DM_PARALLEL
2603 USE module_dm
2604 #endif
2605 ! USE module_date_time
2606 ! USE module_utility
2607
2608
2609 IMPLICIT NONE
2610
2611 ! Arguments
2612 TYPE(domain) :: grid
2613 TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
2614 !INTEGER , INTENT(IN) :: start_step , step , end_step
2615 ! Type (ESMF_Time ) :: start_time, stop_time, CurrTime
2616 ! TYPE(WRFU_TimeInterval) :: time_interval
2617
2618
2619 ! Local data
2620 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
2621 LOGICAL :: emiss_opened
2622 INTEGER :: intime, itime,itime_max, ierr , open_status , fid
2623 REAL :: time, btime, bfrq
2624 REAL, ALLOCATABLE :: dumc0(:,:,:)
2625 CHARACTER (LEN=256) :: message
2626 CHARACTER (LEN=80) :: bdyname
2627
2628 CHARACTER (LEN=9 ),DIMENSION(30) :: ename
2629 INTEGER :: nv,i , j , k, &
2630 ids, ide, jds, jde, kds, kde, &
2631 ims, ime, jms, jme, kms, kme, &
2632 ips, ipe, jps, jpe, kps, kpe
2633
2634
2635 #include <wrf_io_flags.h>
2636
2637 write(message, '(A,I9)') 'call read emissions', intime
2638 call wrf_message( TRIM( message ) )
2639
2640 IF(intime == 0 ) THEN
2641 CALL construct_filename1 ( bdyname , '../../run/wrfem12k_00to12z' , grid%id , 2 )
2642
2643 IF (wrf_dm_on_monitor()) THEN
2644 open (91,file=bdyname,form='unformatted')
2645 ENDIf
2646 write(message, '(A,A)') ' OPENED FILE: ',bdyname
2647 call wrf_message( TRIM( message ) )
2648 ENDIF
2649 IF(intime == 12 ) THEN
2650 CALL construct_filename1 ( bdyname , '../../run/wrfem12k_12to24z' , grid%id , 2 )
2651
2652 IF (wrf_dm_on_monitor()) THEN
2653 open (91,file=bdyname,form='unformatted')
2654 ENDIf
2655 write(message, '(A,A)') ' OPENED FILE: ',bdyname
2656 call wrf_message( TRIM( message ) )
2657 ENDIF
2658 CALL wrf_debug( 100 , 'med_read_bin_chem_emiss: calling emissions' )
2659
2660 CALL get_ijk_from_grid ( grid , &
2661 ids, ide, jds, jde, kds, kde, &
2662 ims, ime, jms, jme, kms, kme, &
2663 ips, ipe, jps, jpe, kps, kpe )
2664
2665 ALLOCATE (dumc0(ids:ide-1,kds:grid%kemit,jds:jde-1))
2666
2667 write(message, '(A,6I6)') ' I am reading emissions, dims: =',ids, ide-1, jds, jde-1, kds, grid%kemit
2668 call wrf_message( TRIM( message ) )
2669
2670 IF(intime == 0 .or. intime == 12) then
2671 read(91)nv
2672 read(91)ename
2673 write(message, '(A,I10)') ' Number of emissions: ',nv
2674 call wrf_message( TRIM( message ) )
2675
2676 ! write(message, '(A,30A10)') ' Array names : ',ename
2677 ! call wrf_message( TRIM( message ) )
2678 ENDIF
2679 read(91)itime
2680 write(message, '(A,I8,A,I8)') ' EMISSIONS INPUT FILE TIME PERIOD (GMT): ',itime-1,' TO ',itime
2681 call wrf_message( TRIM( message ) )
2682
2683 read(91)dumc0
2684 grid%e_so2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2685 read(91)dumc0
2686 grid%e_no(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2687 read(91)dumc0
2688 grid%e_ald(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2689 read(91)dumc0
2690 grid%e_hcho(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2691 read(91)dumc0
2692 grid%e_ora2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2693 read(91)dumc0
2694 grid%e_nh3(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2695 read(91)dumc0
2696 grid%e_hc3(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2697 read(91)dumc0
2698 grid%e_hc5(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2699 read(91)dumc0
2700 grid%e_hc8(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2701 read(91)dumc0
2702 grid%e_eth(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2703 read(91)dumc0
2704 grid%e_co(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2705 read(91)dumc0
2706 grid%e_ol2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2707 read(91)dumc0
2708 grid%e_olt(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2709 read(91)dumc0
2710 grid%e_oli(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2711 read(91)dumc0
2712 grid%e_tol(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2713 read(91)dumc0
2714 grid%e_xyl(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2715 read(91)dumc0
2716 grid%e_ket(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2717 read(91)dumc0
2718 grid%e_csl(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2719 read(91)dumc0
2720 grid%e_iso(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2721 read(91)dumc0
2722 grid%e_pm25i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2723 read(91)dumc0
2724 grid%e_pm25j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2725 read(91)dumc0
2726 grid%e_so4i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2727 read(91)dumc0
2728 grid%e_so4j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2729 read(91)dumc0
2730 grid%e_no3i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2731 read(91)dumc0
2732 grid%e_no3j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2733 read(91)dumc0
2734 grid%e_orgi(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2735 read(91)dumc0
2736 grid%e_orgj(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2737 read(91)dumc0
2738 grid%e_eci(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2739 read(91)dumc0
2740 grid%e_ecj(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2741 read(91)dumc0
2742 grid%e_pm10(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2743
2744 DEALLOCATE ( dumc0 )
2745 RETURN
2746 END SUBROUTINE med_read_bin_chem_emiss
2747
2748 SUBROUTINE med_read_bin_chem_fireemiss ( grid , config_flags ,intime, itime_max)
2749 USE module_bc_time_utilities
2750 #ifdef DM_PARALLEL
2751 USE module_dm
2752 #endif
2753
2754
2755 IMPLICIT NONE
2756
2757 ! Arguments
2758 TYPE(domain) :: grid
2759 TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
2760 !INTEGER , INTENT(IN) :: start_step , step , end_step
2761 ! Type (ESMF_Time ) :: start_time, stop_time, CurrTime
2762 ! TYPE(WRFU_TimeInterval) :: time_interval
2763
2764
2765 ! Local data
2766 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
2767 LOGICAL :: emiss_opened
2768 INTEGER :: intime, itime,itime_max, ierr , open_status , fid
2769 REAL :: time, btime, bfrq
2770 REAL, ALLOCATABLE :: dumc0(:,:)
2771 CHARACTER (LEN=256) :: message
2772 CHARACTER (LEN=80) :: bdyname
2773
2774 CHARACTER (LEN=9 ),DIMENSION(30) :: ename
2775 INTEGER :: nv,i , j , k, &
2776 ids, ide, jds, jde, kds, kde, &
2777 ims, ime, jms, jme, kms, kme, &
2778 ips, ipe, jps, jpe, kps, kpe
2779
2780
2781 #include <wrf_io_flags.h>
2782
2783 write(message, '(A,I9)') 'call read fire emissions', intime
2784 call wrf_message( TRIM( message ) )
2785
2786 IF(intime == 0 ) THEN
2787 CALL construct_filename1 ( bdyname , '../../run/fireemiss' , grid%id , 2 )
2788
2789 IF (wrf_dm_on_monitor()) THEN
2790 open (91,file=bdyname,form='unformatted')
2791 ENDIf
2792 write(message, '(A,A)') ' OPENED FILE: ',bdyname
2793 call wrf_message( TRIM( message ) )
2794 ENDIF
2795 CALL wrf_debug( 100 , 'med_read_bin_chem_fireemiss: calling emissions' )
2796
2797 CALL get_ijk_from_grid ( grid , &
2798 ids, ide, jds, jde, kds, kde, &
2799 ims, ime, jms, jme, kms, kme, &
2800 ips, ipe, jps, jpe, kps, kpe )
2801
2802 ALLOCATE (dumc0(ids:ide-1,jds:jde-1))
2803
2804 write(message, '(A,4I6)') ' I am reading fire emissions, dims: =',ids, ide-1, jds, jde-1
2805 call wrf_message( TRIM( message ) )
2806
2807 IF(intime == 0 .or. intime == 12) then
2808 read(92)nv
2809 read(92)ename
2810 write(message, '(A,I10)') ' Number of emissions: ',nv
2811 call wrf_message( TRIM( message ) )
2812
2813 ENDIF
2814 read(92)itime
2815 write(message, '(A,I8,A,I8)') ' EMISSIONS INPUT FILE TIME PERIOD (GMT): ',itime-1,' TO ',itime
2816 call wrf_message( TRIM( message ) )
2817
2818 read(92)dumc0
2819 ! grid%e_so2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2820 read(92)dumc0
2821 ! grid%e_no(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2822 read(92)dumc0
2823 ! grid%e_ald(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2824 read(92)dumc0
2825 ! grid%e_hcho(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2826 read(92)dumc0
2827 ! grid%e_ora2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2828 read(92)dumc0
2829 ! grid%e_nh3(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2830 read(92)dumc0
2831 ! grid%e_hc3(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2832 read(92)dumc0
2833 ! grid%e_hc5(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2834 read(92)dumc0
2835 ! grid%e_hc8(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2836 read(92)dumc0
2837 ! grid%e_eth(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2838 read(92)dumc0
2839 ! grid%e_co(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2840 read(92)dumc0
2841 ! grid%e_ol2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2842 read(92)dumc0
2843 ! grid%e_olt(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2844 read(92)dumc0
2845 ! grid%e_oli(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2846 read(92)dumc0
2847 ! grid%e_tol(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2848 read(92)dumc0
2849 ! grid%e_xyl(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2850 read(92)dumc0
2851 ! grid%e_ket(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2852 read(92)dumc0
2853 ! grid%e_csl(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2854 read(92)dumc0
2855 ! grid%e_iso(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2856 read(92)dumc0
2857 ! grid%e_pm25i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2858 read(92)dumc0
2859 ! grid%e_pm25j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2860 read(92)dumc0
2861 ! grid%e_so4i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2862 read(92)dumc0
2863 ! grid%e_so4j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2864 read(92)dumc0
2865 ! grid%e_no3i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2866 read(92)dumc0
2867 ! grid%e_no3j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2868 read(92)dumc0
2869 ! grid%e_orgi(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2870 read(92)dumc0
2871 ! grid%e_orgj(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2872 read(92)dumc0
2873 ! grid%e_eci(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2874 read(92)dumc0
2875 ! grid%e_ecj(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2876 read(92)dumc0
2877 ! grid%e_pm10(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2878
2879 DEALLOCATE ( dumc0 )
2880 RETURN
2881 END SUBROUTINE med_read_bin_chem_fireemiss
2882
2883 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2884 END MODULE module_input_chem_data
2885