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