2 !dis Open Source License/Disclaimer, Forecast Systems Laboratory
3 !dis NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
5 !dis This software is distributed under the Open Source Definition,
6 !dis which may be found at http://www.opensource.org/osd.html.
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:
12 !dis - Redistributions of source code must retain this notice, this
13 !dis list of conditions and the following disclaimer.
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.
19 !dis - All modifications to this software must be clearly documented,
20 !dis and are solely the responsibility of the agent making the
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.
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.
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)
48 MODULE module_input_chem_data
52 USE module_driver_constants
53 USE module_state_description
59 USE module_aerosols_sorgam
60 USE module_data_sorgam
62 USE module_get_file_names
67 ! REAL :: grav = 9.8104
68 REAL, PARAMETER :: mwso4 = 96.0576
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
79 PARAMETER( kx=16, kxm1=kx-1, logg=100, lo=34)
81 INTEGER, DIMENSION(logg) :: iref
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
89 DATA so4vaptoaer/.999/
91 CHARACTER (LEN=4), DIMENSION(logg) :: ggnam
93 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
95 ! The idealized profile is converted from the NALROM chemistry model.
97 ! The variable "iref" is the reference index
98 ! "fracref" is the reference fraction correpsonding to iref
100 ! For example: WRF-Chem species number 1 is SO2. iref(1) is 12, and XL(12,K)
101 ! is the profile for SO2.
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).
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.
112 ! The following is a list of long-lived species
118 ! NAMEL( 5)='CH3OOH '
120 ! NAMEL( 7)='ISOPRENE '
122 ! NAMEL( 9)='CH3CHO '
124 ! NAMEL(11)='OTHER ALKA'
126 ! NAMEL(13)='BUTANE '
127 ! NAMEL(14)='ETHENE '
128 ! NAMEL(15)='PROPENE+ '
139 ! NAMEL(26)='PROPANE '
140 ! NAMEL(27)='ACETYLENE '
143 ! NAMEL(30)='NO3+N2O5 '
144 ! NAMEL(31)='HO2NO2 '
145 ! NAMEL(32)='SUM RO2 '
149 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
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/
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./
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' /
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 /
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., &
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., &
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., &
185 #if ( ! NMM_CORE == 0 )
187 DATA ZFA_BDY/ 0., 85., 212., 385., 603., 960., 1430., 2010., &
188 2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
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., &
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/
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/
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/
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/
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/
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/
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/
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/
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/
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/
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/
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/
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/
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/
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/
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/
339 SUBROUTINE input_ext_chem_file (si_grid)
340 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
344 TYPE(domain) :: si_grid
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
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
364 LOGICAL :: input_from_file , need_new_file
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
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
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
381 ! Interface block for routine that passes pointers and needs to know that they
382 ! are receiving pointers.
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.
388 CALL geth_julgmt ( config_flags%julyr , config_flags%julday , config_flags%gmt )
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)
398 write(message,'(A,A)') 'Subroutine input_chem: finding data at date/time: ',extract_date
399 CALL wrf_message ( TRIM(message) )
402 ! And here is an instance of using the information in the NAMELIST.
404 CALL nl_get_debug_level ( 1,debug_level )
405 CALL set_wrf_debug_level ( debug_level )
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.
411 NULLIFY( null_domain )
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...
417 ! if( si_grid%id .EQ. 1) then
418 CALL alloc_and_configure_domain ( domain_id = 1 , &
420 parent = null_domain , &
424 ! CALL alloc_and_configure_domain ( domain_id = 2 , &
425 ! grid = chem_grid , &
426 ! parent = parent_grid , &
432 CALL wrf_debug ( 100 , 'wrfchem:input_chem: set pointer for domain 1' )
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
440 need_new_file = .FALSE.
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) )
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 )
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) )
454 CALL open_r_dataset ( fid, TRIM(inpname) , chgrid, config_flags, "DATASET=INPUT", ierr )
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 )
461 ! How many data time levels in the input file?
465 CALL wrf_debug ( 100, 'subroutine input_chem: find time_loop_max' )
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
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 )
475 write(message,'(6A)') 'Subroutine input_chem: file date/time = ',date_string,' status = ', status_next_var
476 CALL wrf_message ( TRIM(message) )
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 )
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
493 EXIT get_the_right_time
495 CYCLE get_the_right_time
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 )
506 ! For now, the input date and time MUST match
508 ! Put the time check here and set num_steps
510 num_steps(1) = time_loop_max
511 num_steps(2) = time_loop_max+1
512 use_date = date_string
515 EXIT get_the_right_time
518 if( num_steps(2) == time_loop_max ) then
521 END DO get_the_right_time
523 num_steps(2) = MIN(num_steps(2),time_loop_max)
525 ! wgt0 = (ext_time - wgt_time1) / (wgt_time2 - wgt_time1)
528 ! Make sure the right date and time for the chemin data has been found
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 ) )
535 WRITE( wrf_err_message, FMT='(A)' ) 'subroutine input_chem: error finding chemin date/time #2'
536 CALL WRF_ERROR_FATAL ( wrf_err_message )
539 ! There has to be a more elegant way to get to the beginning of the file, but this will do.
541 CALL close_dataset ( fid , config_flags , "DATASET=INPUT" )
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) )
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 )
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.
556 big_time_loop_thingy : DO time_loop = 1 , time_loop_max
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' )
562 IF( time_loop .EQ. num_steps(1) ) THEN
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 )
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) )
576 si_zsigf = (si_grid%em_ph_1 + si_grid%em_phb) / grav
577 ch_zsigf = ( chgrid%em_ph_1 + chgrid%em_phb) / grav
580 #if ( NMM_CORE == 1 )
582 ! Get scalar grid point heights
583 ALLOCATE( pint(ips:ipe,kps:kpe,jps:jpe) )
584 ALLOCATE( pdsl(ips:ipe,jps:jpe) )
586 IF(chgrid%sigma.EQ. 1)THEN
589 pdsl(i,j)=chgrid%nmm_pd(i,j)
595 pdsl(i,j)=chgrid%nmm_res(i,j)*chgrid%nmm_pd(i,j)
600 !!!!! SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
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)
611 CALL wrf_debug (0, message)
613 IF(si_grid%sigma.EQ. 1)THEN
616 pdsl(i,j)=si_grid%nmm_pd(i,j)
622 pdsl(i,j)=si_grid%nmm_res(i,j)*si_grid%nmm_pd(i,j)
626 ! write(message,'(1e15.6,i6)') pdsl(ips,jps), si_grid%sigma
628 !!!!! SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
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
637 si_zsigf(i,k,j)=pint(i,k,j)
642 ! write(message,'(4e15.6)') si_zsigf(1,1:4,1)
643 ! CALL wrf_error_fatal (message)
645 DEALLOCATE( pint); DEALLOCATE( pdsl )
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,:) )
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,:) )
656 ! check size of si_grid vs chgrid
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
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 )
674 ! Fill top level to prevent spurious interpolation results (no extrapolation)
675 chgrid%chem(:,kde,:,:) = chgrid%chem(:,kde-1,:,:)
677 ! Interpolate the chemistry data to the SI grid (holding place for time interpolation)
679 call vinterp_chem(ims, ime, jms, jme, kms, kme, kme, num_chem, ch_zsig, si_zsig, &
680 chgrid%chem, si_grid%chem, .false.)
682 if(wgt0 == 0.) EXIT big_time_loop_thingy
685 IF( time_loop .EQ. num_steps(2) ) THEN
687 ! ! input chemistry sigma levels
688 ! ch_zsigf = ( chgrid%em_ph_1 + chgrid%em_phb) / grav
690 ! ch_zsig(:,k,:) = 0.5 * ( ch_zsigf(:,k,:) + ch_zsigf(:,k+1,:) )
692 ! ch_zsig(:,kde,:) = 0.5 * ( 3. * ch_zsigf(:,kde,:) - ch_zsigf(:,kde-1,:) )
694 ! ! Fill top level to prevent spurious interpolation results (no extrapolation)
695 ! chgrid%chem(:,kde,:,:) = chgrid%chem(:,kde-1,:,:)
697 ! ! Interpolate the chemistry data to the temp chgrid2 structure
699 ! call vinterp_chem(ims, ime, jms, jme, kms, kme, kme, num_chem, ch_zsig, si_zsig, &
700 ! chgrid%chem, chgrid2%chem, .false.)
702 ! ! use linear interpolation in time to get new chem arrays
703 ! si_grid%chem = (1. - wgt0) * si_grid%chem + &
704 ! wgt0 * chgrid2%chem
706 DEALLOCATE( si_zsigf); DEALLOCATE( si_zsig )
707 DEALLOCATE( ch_zsigf); DEALLOCATE( ch_zsig )
709 EXIT big_time_loop_thingy
711 END DO big_time_loop_thingy
713 ! Check for errors in chemin data set
719 si_grid%chem(i,k,j,l) = MAX(si_grid%chem(i,k,j,l),epsilc)
726 ! Close chemin data set
727 CALL close_dataset ( fid , config_flags , "DATASET=INPUT" )
730 CALL domain_destroy( chem_grid )
732 CALL wrf_debug ( 100,' input_ext_chem_data: exit subroutine ')
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)
741 ! Interpolates columns of chemistry data from one set of height surfaces to
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
762 ! Loop over the number of chemical species
763 chem_loop: DO l = 2, nch
765 data_out(:,:,:,l) = -99999.9
770 output_loop: DO k = nz1, nz_out
773 desired_z = z_out(i,k,j)
774 IF (desired_z .LT. z_in(i,1,j)) THEN
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)
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
785 ! Add a check to make sure we are not using the gradient of
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) )
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) )
795 data_out(i,k,j,l) = MAX( data_in(i,1,j,l) + &
796 dvaldz * (desired_z-z_in(i,1,j)), 0.)
798 data_out(i,k,j,l) = data_in(i,1,j,l)
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)
805 IF (extrapolate) THEN
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))
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))
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.)
817 data_out(i,k,j,l) = data_in(i,nz_in,j,l)
821 ! We can trap between two levels and linearly interpolate
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)
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.)
839 #if ( NMM_CORE == 1 )
841 desired_z = z_out(i,k,j)
842 IF (desired_z .GT. z_in(i,1,j)) THEN
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)
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
853 ! Add a check to make sure we are not using the gradient of
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) )
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) )
863 data_out(i,k,j,l) = MAX( data_in(i,1,j,l) + &
864 dvaldz * (desired_z-z_in(i,1,j)), 0.)
866 data_out(i,k,j,l) = data_in(i,1,j,l)
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)
873 IF (extrapolate) THEN
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))
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))
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.)
885 data_out(i,k,j,l) = data_in(i,nz_in,j,l)
889 ! We can trap between two levels and linearly interpolate
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)
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.)
913 END SUBROUTINE vinterp_chem
914 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
915 SUBROUTINE input_chem_profile (si_grid)
919 TYPE(domain) :: si_grid
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
928 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: si_zsigf, si_zsig
930 #if ( ! NMM_CORE == 0 )
931 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: pint
932 REAL, ALLOCATABLE, DIMENSION(:,:) :: pdsl
935 CHARACTER (LEN=80) :: inpname, message
937 write(message,'(A)') 'Subroutine input_chem_profile: '
938 CALL wrf_message ( TRIM(message) )
940 ! And here is an instance of using the information in the NAMELIST.
942 CALL nl_get_debug_level ( 1,debug_level )
943 CALL set_wrf_debug_level ( debug_level )
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 )
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) )
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
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) )
964 write(message,'(A)') 'WRF_NMM_CORE '
965 CALL wrf_message ( message )
967 IF(si_grid%sigma.EQ. 1)THEN
970 pdsl(i,j)=si_grid%nmm_pd(i,j)
976 pdsl(i,j)=si_grid%nmm_res(i,j)*si_grid%nmm_pd(i,j)
984 !!!!! SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
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
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)
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
1004 ! si_zsigf = si_grid%z
1007 ! si_zsigf = (si_grid%em_ph_1 + si_grid%em_phb) / grav
1010 si_zsig(:,k,:) = 0.5 * ( si_zsigf(:,k,:) + si_zsigf(:,k+1,:) )
1012 si_zsig(:,kde,:) = 0.5 * ( 3. * si_zsigf(:,kde,:) - si_zsigf(:,kde-1,:) )
1014 ! An alternative ozone profile option for initialization
1015 if( si_grid%gas_ic_opt == GAS_IC_PNNL ) &
1016 call gasprofile_init_pnnl
1018 ! Determine the index of the last gas species
1019 numgas = get_last_gas(si_grid%chem_opt)
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
1029 CALL make_chem_profile (ims, ime, jms, jme, kms, kme, num_chem, numgas, &
1030 si_zsig, si_grid%chem)
1033 CALL wrf_debug ( 100,' input_chem_profile: exit subroutine ')
1035 DEALLOCATE( si_zsigf ); DEALLOCATE( si_zsig )
1036 #if ( ! NMM_CORE == 0 )
1037 DEALLOCATE( pdsl ); DEALLOCATE( pint )
1041 END SUBROUTINE input_chem_profile
1042 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1043 SUBROUTINE make_chem_profile ( nx1, nx2, ny1, ny2, nz1, nz2, nch, numgas, &
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
1053 CHARACTER (LEN=80) :: message
1054 INTEGER :: i, j, k, l, is
1056 REAL, DIMENSION(nx1:nx2,nz1:kx ,ny1:ny2,lo+1):: chprof
1057 REAL, DIMENSION(nx1:nx2,nz1:kx ,ny1:ny2) :: zprof
1059 REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2,nch) :: chem
1060 REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2,lo ) :: stor
1061 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1063 ! Check the number of species
1065 if( nch .NE. num_chem) then
1066 message = ' Input_chem_profile: wrong number of chemical species'
1067 ! CALL WRF_ERROR_FATAL ( message )
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.
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))
1082 ! return xl to previous value for next time...
1083 ! 34 chemicals (lo), 16 vertical levels (kx)
1085 ! xl(i,1:kx)=xl(i,1:kx)*dens(1:kx)
1088 ! Change number concentrations to mixing ratios for short-lived NALROM species
1090 chprof(:,k,:,lo-5:lo+1) = chprof(:,k,:,lo-5:lo+1)/dens(k)
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.)
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)
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.
1110 chem(i,k,j,l)=fracref(l-1)*stor(i,k,j,is)*1.E6
1117 END SUBROUTINE make_chem_profile
1118 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1120 ! this is a kludge routine as of now....
1122 SUBROUTINE bdy_chem_value_sorgam (chem, z, nch, config_flags, &
1124 USE module_data_sorgam
1125 USE module_aerosols_sorgam
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
1136 REAL, DIMENSION(lo+1,1:kx):: cprof ! chemical profile, diff. index order
1138 REAL, DIMENSION(1:kx):: zprof
1139 REAL, DIMENSION(lo ) :: stor
1142 real :: chemsulf_radm,chem_so4aj,chem_so4ai
1145 !between gas and aerosol phase
1147 !factor for splitting initial conc. of SO4
1148 !3rd moment i-mode [3rd moment/m^3]
1150 !3rd MOMENT j-mode [3rd moment/m^3]
1157 ! method for bc calculation is determined by aer_bc_opt
1159 if (config_flags%aer_bc_opt == AER_BC_PNNL) then
1160 call sorgam_set_aer_bc_pnnl( chem, z, nch )
1162 else if (config_flags%aer_bc_opt == AER_BC_DEFAULT) then
1165 call wrf_error_fatal( &
1166 "bdy_chem_value_sorgam -- unable to parse aer_bc_opt" )
1169 ! do default calculation of sorgam aerosol bc values
1171 ! tempfac=(t+t0)*((p+pb)/p1000mb)**rcp
1172 ! convfac=(p+pb)/rgasuniv/tempfac
1174 !--- units for advection....
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
1183 ! recalculate sulf profile for aerosols
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
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)
1196 zprof(k) = 0.5*(zfa_bdy(k)+zfa_bdy(k+1))
1198 cprof(l+1,k) = xl(l,kx+1-k)
1200 ! Fix number concentrations to mixing ratios for short-lived NALROM species
1202 cprof(l+1,k) = xl(l,kx+1-k)/dens(kx+1-k)
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)
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)
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.)
1228 ! Here is where the chemistry value is constructed
1229 chemsulf_radm = fracref(p_sulf-1)*stor( iref(p_sulf-1) )*1.E6
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)
1241 ! compute values for aerosol input data
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)
1248 l2sginin = xxlsgn**2
1249 l2sginia = xxlsga**2
1250 l2sginic = xxlsgc**2
1252 en1 = exp(0.125*l2sginin)
1253 ea1 = exp(0.125*l2sginia)
1254 ec1 = exp(0.125*l2sginic)
1271 esn12 = esn04*esn04*esn04
1272 esa12 = esa04*esa04*esa04
1273 esc12 = esc04*esc04*esc04
1304 ! Units are something like number concentration
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
1312 END SUBROUTINE bdy_chem_value_sorgam
1314 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1315 SUBROUTINE bdy_chem_value_tracer ( chem, nch )
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
1325 REAL, intent(OUT) :: chem
1326 INTEGER, intent(IN) :: nch ! index number of chemical species
1327 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1329 if( nch .ne. p_co ) then
1331 else if( nch == p_co ) then
1337 END SUBROUTINE bdy_chem_value_tracer
1338 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1340 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1341 SUBROUTINE bdy_chem_value_racm ( chem, z, nch, numgas,p_co2 )
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
1350 INTEGER :: i, k, irefcur
1352 REAL, DIMENSION(kx):: cprof ! chemical profile, diff. index order
1354 REAL, DIMENSION(1:kx):: zprof
1358 CHARACTER (LEN=80) :: message
1359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1361 ! Check the number of species
1362 ! if((nch-1).gt.logg)return
1363 if (nch.eq.p_co2)then
1367 if (nch.eq.p_co2+1)then
1371 if (nch.ge.p_co2+2)return
1374 ! if( nch .GT. logg+1) then
1375 if( nch .GT. numgas) then
1376 message = ' Input_chem_profile: wrong number of chemical species'
1378 ! CALL WRF_ERROR_FATAL ( message )
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)
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)
1391 cprof(k) = xl(irefcur,kx+1-k)/dens(kx+1-k)
1395 ! Interpolate temp 3D chemical profile array to WRF field
1396 IF (z .LT. zprof(1)) THEN
1398 ELSE IF (z .GT. zprof(kx)) THEN
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
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.)
1417 ! Here is where the chemistry value is constructed
1418 chem = fracref(nch-1)*stor*1.E6
1420 ! special code for sulfate/h2so4
1421 if(nch.eq.p_sulf.and.p_nu0.gt.1)then
1422 chem=chem*(1.-so4vaptoaer)
1426 END SUBROUTINE bdy_chem_value_racm
1427 SUBROUTINE bdy_chem_value ( chem, z, nch, numgas )
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
1436 INTEGER :: i, k, irefcur
1438 REAL, DIMENSION(kx):: cprof ! chemical profile, diff. index order
1440 REAL, DIMENSION(1:kx):: zprof
1444 CHARACTER (LEN=80) :: message
1445 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1447 ! Check the number of species
1448 ! if((nch-1).gt.logg)return
1450 ! if( nch .GT. logg+1) then
1451 if( nch .GT. numgas) then
1452 message = ' Input_chem_profile: wrong number of chemical species'
1454 ! CALL WRF_ERROR_FATAL ( message )
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)
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)
1467 cprof(k) = xl(irefcur,kx+1-k)/dens(kx+1-k)
1471 ! Interpolate temp 3D chemical profile array to WRF field
1472 IF (z .LT. zprof(1)) THEN
1474 ELSE IF (z .GT. zprof(kx)) THEN
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
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.)
1493 ! Here is where the chemistry value is constructed
1494 chem = fracref(nch-1)*stor*1.E6
1496 ! special code for sulfate/h2so4
1497 if(nch.eq.p_sulf.and.p_nu0.gt.1)then
1498 chem=chem*(1.-so4vaptoaer)
1502 END SUBROUTINE bdy_chem_value
1503 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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, &
1514 u, v, config_flags, alt, &
1515 t,pb,p,t0,p1000mb,rcp,ph,phb,g, &
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 )
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.
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
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 ) , &
1548 real, INTENT (IN) :: g,rcp,t0,p1000mb
1549 TYPE( grid_config_rec_type ) config_flags
1551 INTEGER :: i, j, k, numgas
1552 INTEGER :: ibs, ibe, jbs, jbe, itf, jtf, ktf
1553 INTEGER :: i_inner, j_inner
1555 integer :: itestbc, i_bdy_method
1556 real tempfac,convfac
1558 logical :: have_bcs_chem
1560 chem_bv_def = conmin
1561 numgas = get_last_gas(config_flags%chem_opt)
1563 if(p_nu0.gt.1)itestbc=1
1566 itf = min(ite,ide-1)
1569 jtf = min(jte,jde-1)
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
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)
1585 if ((ic .ge. p_so2) .and. (ic .le. p_ho2)) then
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
1596 else if ((ic .ge. p_so4aj) .and. (ic .le. p_corn)) then
1598 else if ((ic .ge. p_hcl) .and. (ic .le. p_isopo2)) then
1600 else if ((ic .ge. p_dms) .and. (ic .le. p_mtf)) then
1602 else if ((ic .ge. p_so4_a01) .and. (ic .le. p_num_a01)) then
1604 else if ((ic .ge. p_so4_a02) .and. (ic .le. p_num_a02)) then
1606 else if ((ic .ge. p_so4_a03) .and. (ic .le. p_num_a03)) then
1608 else if ((ic .ge. p_so4_a04) .and. (ic .le. p_num_a04)) then
1610 else if ((ic .ge. p_so4_a05) .and. (ic .le. p_num_a05)) then
1612 else if ((ic .ge. p_so4_a06) .and. (ic .le. p_num_a06)) then
1614 else if ((ic .ge. p_so4_a07) .and. (ic .le. p_num_a07)) then
1616 else if ((ic .ge. p_so4_a08) .and. (ic .le. p_num_a08)) then
1618 else if (config_flags%chem_opt == CHEM_TRACER) then
1621 if (have_bcs_chem) i_bdy_method =6
1622 if (ic .lt. param_first_scalar) i_bdy_method = 0
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
1636 ! print 90010, '_bdy_NONE** for ic=', ic, i_bdy_method
1638 !90010 format( a, 2(1x,i5) )
1639 !90020 format( a, 1p, 2e12.2 )
1640 !----------------------------------------------------------------------
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
1646 IF (jts - jbs .lt. spec_zone) THEN
1648 DO j = jts, min(jtf,jbs+spec_zone-1)
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)
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
1688 chem(i,k,j) = chem_bv_def
1695 IF (jbe - jtf .lt. spec_zone) THEN
1697 DO j = max(jts,jbe-spec_zone+1), jtf
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)
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)
1729 chem(i,k,j) = chem_bv_def
1737 IF (its - ibs .lt. spec_zone) THEN
1739 DO i = its, min(itf,ibs+spec_zone-1)
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)
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)
1771 chem(i,k,j) = chem_bv_def
1779 IF (ibe - itf .lt. spec_zone) THEN
1781 DO i = max(its,ibe-spec_zone+1), itf
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)
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)
1813 chem(i,k,j) = chem_bv_def
1821 END SUBROUTINE flow_dep_bdy_chem
1823 SUBROUTINE flow_dep_bdy_chem ( chem, chem_b,chem_bt,dt, &
1825 ijds, ijde,have_bcs_chem, &
1826 u, v, config_flags, alt, &
1827 t,pb,p,t0,p1000mb,rcp,ph,phb,g, &
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 )
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.
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
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 ) , &
1861 real, INTENT (IN) :: g,rcp,t0,p1000mb
1862 TYPE( grid_config_rec_type ) config_flags
1864 INTEGER :: i, j, k, numgas
1865 INTEGER :: ibs, ibe, jbs, jbe, itf, jtf, ktf
1866 INTEGER :: i_inner, j_inner
1868 integer :: itestbc, i_bdy_method
1869 real tempfac,convfac
1871 logical :: have_bcs_chem
1873 chem_bv_def = conmin
1874 numgas = get_last_gas(config_flags%chem_opt)
1876 if(p_nu0.gt.1)itestbc=1
1879 itf = min(ite,ide-1)
1882 jtf = min(jte,jde-1)
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
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)
1898 if ((ic .ge. p_so2) .and. (ic .le. p_ho2)) then
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
1909 else if ((ic .ge. p_so4aj) .and. (ic .le. p_corn)) then
1911 else if ((ic .ge. p_hcl) .and. (ic .le. p_isopo2)) then
1913 else if ((ic .ge. p_dms) .and. (ic .le. p_mtf)) then
1915 else if ((ic .ge. p_so4_a01) .and. (ic .le. p_num_a01)) then
1917 else if ((ic .ge. p_so4_a02) .and. (ic .le. p_num_a02)) then
1919 else if ((ic .ge. p_so4_a03) .and. (ic .le. p_num_a03)) then
1921 else if ((ic .ge. p_so4_a04) .and. (ic .le. p_num_a04)) then
1923 else if ((ic .ge. p_so4_a05) .and. (ic .le. p_num_a05)) then
1925 else if ((ic .ge. p_so4_a06) .and. (ic .le. p_num_a06)) then
1927 else if ((ic .ge. p_so4_a07) .and. (ic .le. p_num_a07)) then
1929 else if ((ic .ge. p_so4_a08) .and. (ic .le. p_num_a08)) then
1931 else if (config_flags%chem_opt == CHEM_TRACER) then
1934 if (have_bcs_chem) i_bdy_method =6
1935 if (ic .lt. param_first_scalar) i_bdy_method = 0
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
1949 ! print 90010, '_bdy_NONE** for ic=', ic, i_bdy_method
1951 !90010 format( a, 2(1x,i5) )
1952 !90020 format( a, 1p, 2e12.2 )
1953 !----------------------------------------------------------------------
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
1959 IF (jts - jbs .lt. spec_zone) THEN
1961 DO j = jts, min(jtf,jbs+spec_zone-1)
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)
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
2001 chem(i,k,j) = chem_bv_def
2008 IF (jbe - jtf .lt. spec_zone) THEN
2010 DO j = max(jts,jbe-spec_zone+1), jtf
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)
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)
2042 chem(i,k,j) = chem_bv_def
2050 IF (its - ibs .lt. spec_zone) THEN
2052 DO i = its, min(itf,ibs+spec_zone-1)
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)
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)
2084 chem(i,k,j) = chem_bv_def
2092 IF (ibe - itf .lt. spec_zone) THEN
2094 DO i = max(its,ibe-spec_zone+1), itf
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)
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)
2126 chem(i,k,j) = chem_bv_def
2134 END SUBROUTINE flow_dep_bdy_chem
2136 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2137 SUBROUTINE bdy_chem_value_gcm ( chem, chem_b, chem_bt, dt,ic)
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
2148 CHARACTER (LEN=80) :: message
2149 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2151 ! if( nch .GT. numgas) then
2152 ! message = ' Input_chem_profile: wrong number of chemical species'
2154 ! CALL WRF_ERROR_FATAL ( message )
2158 !print*,'before',chem,chem_bt ,dt, ic
2160 chem=max(epsilc,chem_b + chem_bt * dt)
2161 !print*,'after',chem
2163 END SUBROUTINE bdy_chem_value_gcm
2164 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2166 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2167 SUBROUTINE cv_mmdd_jday ( YY, MM, DD, JDAY)
2169 ! Subroutine to compute the julian day given the month and day
2172 INTEGER, INTENT(IN ) :: YY, MM, DD
2173 INTEGER, INTENT(OUT) :: JDAY
2175 INTEGER, DIMENSION(12) :: imon, imon_a
2178 DATA imon_a /0,31,59,90,120,151,181,212,243,273,304,334/
2180 !..... Check for leap year.
2185 if(YY .eq. (YY/4)*4) then
2187 imon(i) = imon(i) + 1
2191 !..... Convert month, day to julian day.
2193 jday = imon(mm) + dd
2196 END SUBROUTINE cv_mmdd_jday
2198 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2200 integer FUNCTION get_last_gas(chem_opt)
2202 integer, intent(in) :: chem_opt
2204 ! Determine the index of the last gas species, which depends
2205 ! upon the gas mechanism.
2207 select case (chem_opt)
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
2214 get_last_gas = p_mtf
2216 case (CBMZ_BB,CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN,CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ)
2217 get_last_gas = p_isopo2
2223 call wrf_error_fatal("get_last_gas: could not decipher chem_opt value")
2227 END FUNCTION get_last_gas
2228 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2230 !****************************************************************
2232 ! SUBROUTINE TO SET AEROSOL BC VALUES USING THE *
2233 ! aer_bc_opt == aer_bc_pnnl OPTION. *
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 *
2240 ! CALLS THE FOLLOWING SUBROUTINES: NONE *
2242 ! CALLED BY: bdy_chem_value_sorgam *
2244 !****************************************************************
2245 SUBROUTINE sorgam_set_aer_bc_pnnl( chem, z, nch )
2246 USE module_data_sorgam
2249 INTEGER,INTENT(IN ) :: nch
2250 real,intent(in ) :: z
2251 REAL,INTENT(INOUT ) :: chem
2254 m3acc, m3cor, m3nuc, &
2255 bv_so4ai, bv_so4aj, &
2256 bv_nh4ai, bv_nh4aj, &
2257 bv_no3ai, bv_no3aj, &
2260 bv_orgpai,bv_orgpaj, &
2261 bv_antha, bv_seas, bv_soila
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 !!$! --------- ----------
2270 !!$! 2000<z<3000 linear transition zone to 0.5
2271 !!$! 3000<z<5000 linear transision zone to 0.25
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
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.)
2289 ! Updated aerosol profile multiplier 1-Apr-2005:
2290 ! Height(m) Multiplier
2291 ! --------- ----------
2293 ! 2000<z<3000 linear transition zone to 0.25
2294 ! 3000<z<5000 linear transision zone to 0.125
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
2304 .and. z <= 3000. ) then
2305 mult = 1.0 - 0.00075*(z-2000.)
2307 .and. z <= 5000. ) then
2308 mult = 0.25 - 4.166666667e-5*(z-3000.)
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
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 + &
2335 orgfac*8.0*conmin + orgfac*bv_orgpai + &
2336 anthfac*bv_p25i + anthfac*bv_eci
2338 !... j-mode (note that the 8 SOA species have bv=conmin)
2339 m3acc = so4fac*bv_so4aj + nh4fac*bv_nh4aj + &
2341 orgfac*8.0*conmin + orgfac*bv_orgpaj + &
2342 anthfac*bv_p25j + anthfac*bv_ecj
2345 m3cor = soilfac*bv_soila + seasfac*bv_seas + &
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
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
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
2383 if( nch == p_antha ) chem = bv_antha
2384 if( nch == p_soila ) chem = bv_soila
2385 if( nch == p_seas ) chem = bv_seas
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)
2391 END SUBROUTINE sorgam_set_aer_bc_pnnl
2392 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2394 !****************************************************************
2396 ! SUBROUTINE TO OVERWRITE THE PREDEFINED OZONE PROFILE *
2397 ! WHEN gas_ic_opt == gas_ic_pnnl *
2398 ! OR gas_bc_opt == gas_bc_pnnl *
2400 ! wig, 21-Apr-2004 *
2401 ! rce 25-apr-2004 - changed name to *
2402 ! "gasprofile_init_pnnl" *
2404 ! CALLS THE FOLLOWING SUBROUTINES: NONE *
2406 ! CALLED BY: chem_init *
2407 ! input_chem_profile *
2409 !****************************************************************
2410 SUBROUTINE gasprofile_init_pnnl
2411 use module_data_sorgam,only: conmin
2415 call wrf_debug ( 500 , 'wrfchem:gasprofile_init_pnnl' )
2416 ! print*,'gasprofile_init_pnnl redefining o3 and sulf profiles.'
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/
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
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
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
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
2459 if( p_sulf > 1 ) then
2460 xl(iref(p_sulf-1),:) = conmin
2463 end SUBROUTINE gasprofile_init_pnnl
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, &
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, &
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, &
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 ), &
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, &
2514 print*,"itimestep =",itimestep
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)
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, &
2528 print*,"e_org i,j:",e_orgi(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, &
2556 print*,"e_no3 f,c:",e_no3j(i,k,j)*conva, &
2558 print*,"e_orgc:",e_orgc(i,k,j)*conva
2559 print*,"e_ecc:",e_ecc(i,k,j)*conva
2562 print*,"dz8w=0 so cannot show adjusted emissions"
2564 print*,"CHEM_DBG PRINT (PPM or ug/m^3) AT (i,k,j):",i,k,j
2566 print*,n,chem(i,k,j,n)
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)
2596 end subroutine chem_dbg
2599 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2601 SUBROUTINE med_read_bin_chem_emiss ( grid , config_flags ,intime, itime_max)
2604 ! USE module_io_domain
2607 ! USE module_configure
2608 USE module_bc_time_utilities
2612 ! USE module_date_time
2613 ! USE module_utility
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
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
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
2642 #include <wrf_io_flags.h>
2644 write(message, '(A,I9)') 'call read emissions', intime
2645 call wrf_message( TRIM( message ) )
2647 IF(intime == 0 ) THEN
2648 CALL construct_filename1 ( bdyname , '../../run/wrfem12k_00to12z' , grid%id , 2 )
2650 IF (wrf_dm_on_monitor()) THEN
2651 open (91,file=bdyname,form='unformatted')
2653 write(message, '(A,A)') ' OPENED FILE: ',bdyname
2654 call wrf_message( TRIM( message ) )
2656 IF(intime == 12 ) THEN
2657 CALL construct_filename1 ( bdyname , '../../run/wrfem12k_12to24z' , grid%id , 2 )
2659 IF (wrf_dm_on_monitor()) THEN
2660 open (91,file=bdyname,form='unformatted')
2662 write(message, '(A,A)') ' OPENED FILE: ',bdyname
2663 call wrf_message( TRIM( message ) )
2665 CALL wrf_debug( 100 , 'med_read_bin_chem_emiss: calling emissions' )
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 )
2672 ALLOCATE (dumc0(ids:ide-1,kds:grid%kemit,jds:jde-1))
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 ) )
2677 IF(intime == 0 .or. intime == 12) then
2680 write(message, '(A,I10)') ' Number of emissions: ',nv
2681 call wrf_message( TRIM( message ) )
2683 ! write(message, '(A,30A10)') ' Array names : ',ename
2684 ! call wrf_message( TRIM( message ) )
2687 write(message, '(A,I8,A,I8)') ' EMISSIONS INPUT FILE TIME PERIOD (GMT): ',itime-1,' TO ',itime
2688 call wrf_message( TRIM( message ) )
2691 grid%e_so2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2693 grid%e_no(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2695 grid%e_ald(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2697 grid%e_hcho(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2699 grid%e_ora2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2701 grid%e_nh3(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2703 grid%e_hc3(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2705 grid%e_hc5(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2707 grid%e_hc8(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2709 grid%e_eth(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2711 grid%e_co(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2713 grid%e_ol2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2715 grid%e_olt(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2717 grid%e_oli(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2719 grid%e_tol(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2721 grid%e_xyl(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2723 grid%e_ket(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2725 grid%e_csl(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2727 grid%e_iso(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2729 grid%e_pm25i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2731 grid%e_pm25j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2733 grid%e_so4i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2735 grid%e_so4j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2737 grid%e_no3i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2739 grid%e_no3j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2741 grid%e_orgi(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2743 grid%e_orgj(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2745 grid%e_eci(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2747 grid%e_ecj(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2749 grid%e_pm10(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
2751 DEALLOCATE ( dumc0 )
2753 END SUBROUTINE med_read_bin_chem_emiss
2755 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2756 END MODULE module_input_chem_data