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