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