Line data Source code
1 : module aircraft_emit
2 : !-----------------------------------------------------------------------
3 : !
4 : ! Purpose:
5 : ! Manages reading and interpolation of aircraft aerosols
6 : !
7 : ! Authors: Chih-Chieh (Jack) Chen and Cheryl Craig -- February 2010
8 : !
9 : !-----------------------------------------------------------------------
10 : use perf_mod, only : t_startf, t_stopf
11 :
12 : use shr_kind_mod, only : r8 => shr_kind_r8
13 : use cam_abortutils, only : endrun
14 : use spmd_utils, only : masterproc
15 : use tracer_data, only : trfld, trfile
16 : use cam_logfile, only : iulog
17 :
18 : implicit none
19 : private
20 : save
21 :
22 : public :: aircraft_emit_init
23 : public :: aircraft_emit_adv
24 : public :: aircraft_emit_register
25 : public :: aircraft_emit_readnl
26 : public :: get_aircraft
27 :
28 : type :: forcing_air
29 : real(r8) :: mw
30 : character(len=256) :: filelist
31 : character(len=256) :: filename
32 : real(r8), pointer :: times(:)
33 : real(r8), pointer :: levi(:)
34 : character(len=11) :: species
35 : character(len=8) :: units
36 : integer :: nsectors
37 : character(len=32),pointer :: sectors(:)
38 : type(trfld),pointer :: fields(:)
39 : type(trfile) :: file
40 : end type forcing_air
41 :
42 : type(forcing_air),allocatable :: forcings_air(:)
43 :
44 : integer, parameter :: N_AERO = 13
45 : character(len=13) :: aero_names(N_AERO) = (/'ac_SLANT_DIST','ac_TRACK_DIST','ac_HC ','ac_NOX ','ac_PMNV ',&
46 : 'ac_PMSO ','ac_PMFO ','ac_FUELBURN ','ac_CO2 ','ac_H2O ',&
47 : 'ac_SOX ','ac_CO ','ac_BC '/)
48 :
49 : real(r8), parameter :: molmass(N_AERO) = 1._r8
50 :
51 : logical :: advective_tracer(N_AERO) = .false.
52 : character(len=3) :: mixtype(N_AERO) = 'wet'
53 :
54 : real(r8) :: cptmp = 666.0_r8
55 : real(r8) :: qmin = 0.0_r8
56 : logical :: cam_outfld = .false.
57 :
58 : integer :: index_map(N_AERO)
59 : character(len=256) :: air_specifier(N_AERO)=''
60 : character(len=256) :: air_datapath=''
61 : character(len=24) :: air_type = 'CYCLICAL_LIST' ! 'CYCLICAL_LIST'
62 :
63 : logical :: rmv_file = .false.
64 :
65 : integer :: number_flds
66 :
67 : integer :: aircraft_cnt = 0
68 : character(len=16) :: spc_name_list(N_AERO)
69 : character(len=256) :: spc_flist(N_AERO),spc_fname(N_AERO)
70 : logical :: dist(N_AERO)
71 :
72 : contains
73 :
74 1489176 : subroutine get_aircraft(cnt, spc_name_list_out)
75 : integer, intent(out) :: cnt
76 : character(len=16), optional, intent(out) :: spc_name_list_out(N_AERO)
77 : integer :: i
78 :
79 20848464 : spc_name_list_out = ''
80 :
81 1489176 : cnt = aircraft_cnt
82 1489176 : if( cnt>0 ) then
83 0 : do i=1,cnt
84 0 : spc_name_list_out(i) = spc_name_list(i)
85 : end do
86 : end if
87 :
88 1489176 : end subroutine get_aircraft
89 :
90 1536 : subroutine aircraft_emit_register()
91 :
92 : !------------------------------------------------------------------
93 : ! **** Add the aircraft aerosol data to the physics buffer ****
94 : !------------------------------------------------------------------
95 : use ppgrid, only: pver,pcols
96 : use physics_buffer, only : pbuf_add_field, dtype_r8
97 : use tracer_data, only: incr_filename
98 : use constituents, only: cnst_add
99 :
100 : integer :: i,idx, mm, ind, n
101 : character(len=16) :: spc_name
102 : character(len=256) :: filelist, curr_filename
103 : character(len=128) :: long_name
104 : logical :: has_fixed_ubc=.false.
105 : logical :: read_iv=.false.
106 :
107 : !------------------------------------------------------------------
108 : ! Return if air_specifier is blank (no aircraft data to process)
109 : !------------------------------------------------------------------
110 1536 : dist(:) = .false.
111 1536 : aircraft_cnt = 0
112 1536 : if (air_specifier(1) == "") return
113 :
114 : ! count aircraft emission species used in the simulation
115 0 : count_emis: do n=1,N_AERO
116 :
117 0 : if( len_trim(air_specifier(n) ) == 0 ) then
118 : exit count_emis
119 : endif
120 :
121 0 : i = scan(air_specifier(n),'->')
122 0 : spc_name = trim(adjustl(air_specifier(n)(:i-1)))
123 0 : filelist = trim(adjustl(air_specifier(n)(i+2:)))
124 :
125 0 : mm = get_aircraft_ndx(spc_name)
126 0 : if( mm < 1 ) then
127 0 : call endrun('aircraft_emit_register: '//trim(spc_name)//' is not in the aircraft emission dataset')
128 : endif
129 :
130 0 : if (trim(spc_name) == 'ac_SLANT_DIST'.or. trim(spc_name) == 'ac_TRACK_DIST') dist(n) = .true.
131 :
132 0 : aircraft_cnt = aircraft_cnt + 1
133 0 : call pbuf_add_field(aero_names(mm),'physpkg',dtype_r8,(/pcols,pver/),idx)
134 :
135 0 : spc_flist(aircraft_cnt) = filelist
136 0 : spc_name_list(aircraft_cnt) = spc_name
137 0 : index_map(aircraft_cnt) = mm
138 :
139 0 : curr_filename=''
140 : spc_fname(aircraft_cnt) = incr_filename( curr_filename, filenames_list=spc_flist(aircraft_cnt), &
141 0 : datapath=air_datapath)
142 :
143 0 : if( advective_tracer(mm) ) then
144 0 : long_name = 'aircraft_'//trim(spc_name)
145 : call cnst_add(aero_names(mm),molmass(mm),cptmp,qmin,ind,longname=long_name,readiv=read_iv, &
146 0 : mixtype=mixtype(mm),cam_outfld=cam_outfld,fixed_ubc=has_fixed_ubc)
147 : endif
148 :
149 : enddo count_emis
150 : ! count aircraft emission species used in the simulation
151 :
152 1536 : endsubroutine aircraft_emit_register
153 :
154 1536 : subroutine aircraft_emit_init()
155 : !-------------------------------------------------------------------
156 : ! **** Initialize the aircraft aerosol data handling ****
157 : !-------------------------------------------------------------------
158 1536 : use cam_history, only: addfld, add_default
159 : use tracer_data, only: trcdata_init
160 : use phys_control, only: phys_getopts
161 :
162 : implicit none
163 :
164 : character(len=16) :: spc_name
165 :
166 : integer :: astat, m
167 :
168 : logical :: history_chemistry
169 :
170 1536 : call phys_getopts(history_chemistry_out=history_chemistry)
171 :
172 : !------------------------------------------------------------------
173 : ! Return if aircraft_cnt is zero (no aircraft data to process)
174 : !------------------------------------------------------------------
175 1536 : if (aircraft_cnt == 0 ) return
176 :
177 0 : if (masterproc) write(iulog,*) ' '
178 :
179 : !-----------------------------------------------------------------------
180 : ! allocate forcings type array
181 : !-----------------------------------------------------------------------
182 0 : allocate( forcings_air(aircraft_cnt), stat=astat )
183 0 : if( astat/= 0 ) then
184 0 : write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air array; error = ',astat
185 0 : call endrun
186 : end if
187 :
188 : !-----------------------------------------------------------------------
189 : ! setup the forcings_air type array
190 : !-----------------------------------------------------------------------
191 0 : species_loop : do m = 1,aircraft_cnt
192 :
193 0 : allocate( forcings_air(m)%sectors(1), stat=astat )
194 0 : if( astat/= 0 ) then
195 0 : write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air%sectors array; error = ',astat
196 0 : call endrun
197 : end if
198 :
199 0 : allocate( forcings_air(m)%fields(1), stat=astat )
200 0 : if( astat/= 0 ) then
201 0 : write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air%fields array; error = ',astat
202 0 : call endrun
203 : end if
204 :
205 0 : spc_name = spc_name_list(m)
206 : !-----------------------------------------------------------------------
207 : ! default settings
208 : !-----------------------------------------------------------------------
209 0 : forcings_air(m)%file%stepTime = .true. ! Aircraft data is not to be interpolated in time
210 0 : forcings_air(m)%file%cyclical_list = .true. ! Aircraft data cycles over the filename list
211 0 : forcings_air(m)%file%weight_by_lat = .true. ! Aircraft data - interpolated with latitude weighting
212 0 : forcings_air(m)%file%conserve_column = .true. ! Aircraft data - vertically interpolated to conserve the total column
213 0 : forcings_air(m)%file%dist = dist(m)
214 0 : forcings_air(m)%species = spc_name
215 0 : forcings_air(m)%sectors = spc_name ! Only one species per file for aircraft data
216 0 : forcings_air(m)%nsectors = 1
217 0 : forcings_air(m)%filelist = spc_flist(m)
218 : ! forcings_air(m)%file%curr_filename = spc_fname(m)
219 0 : forcings_air(m)%filename = spc_fname(m)
220 : end do species_loop
221 :
222 0 : if (masterproc) then
223 : !-----------------------------------------------------------------------
224 : ! diagnostics
225 : !-----------------------------------------------------------------------
226 0 : write(iulog,*) ' '
227 0 : write(iulog,*) 'aircraft_emit_init: diagnostics'
228 0 : write(iulog,*) ' '
229 0 : write(iulog,*) 'aircraft_emit timing specs'
230 0 : write(iulog,*) 'type = ',air_type
231 0 : write(iulog,*) ' '
232 0 : write(iulog,*) 'there are ',aircraft_cnt,' species of aircraft emission'
233 0 : do m = 1,aircraft_cnt
234 0 : write(iulog,*) ' '
235 0 : write(iulog,*) 'forcing type ',m
236 0 : write(iulog,*) 'species = ',trim(forcings_air(m)%species)
237 0 : write(iulog,*) 'filelist= ',trim(forcings_air(m)%filelist)
238 : end do
239 0 : write(iulog,*) ' '
240 : endif
241 :
242 : !------------------------------------------------------------------
243 : ! Initialize the aircraft file processing
244 : !------------------------------------------------------------------
245 0 : do m=1,aircraft_cnt
246 :
247 0 : allocate (forcings_air(m)%file%in_pbuf(size(forcings_air(m)%sectors)))
248 0 : forcings_air(m)%file%in_pbuf(:) = .true.
249 0 : call trcdata_init( forcings_air(m)%sectors, forcings_air(m)%filename, forcings_air(m)%filelist, air_datapath, &
250 0 : forcings_air(m)%fields, forcings_air(m)%file, rmv_file, 0, 0, 0, air_type)
251 :
252 :
253 0 : number_flds = 0
254 0 : if (associated(forcings_air(m)%fields)) number_flds = size( forcings_air(m)%fields )
255 :
256 0 : if( number_flds < 1 ) then
257 0 : if ( masterproc ) then
258 0 : write(iulog,*) 'There are no aircraft aerosols'
259 0 : write(iulog,*) ' '
260 0 : call endrun
261 : endif
262 : end if
263 :
264 0 : spc_name = spc_name_list(m)
265 0 : call addfld( trim(spc_name), (/ 'lev' /), 'A', forcings_air(m)%fields(1)%units, &
266 0 : 'aircraft emission '//trim(spc_name) )
267 0 : if (history_chemistry) then
268 0 : call add_default( trim(spc_name), 1, ' ' )
269 : end if
270 : end do
271 :
272 :
273 3072 : end subroutine aircraft_emit_init
274 :
275 :
276 :
277 741888 : subroutine aircraft_emit_adv( state, pbuf2d)
278 : !-------------------------------------------------------------------
279 : ! **** Advance to the next aircraft data ****
280 : !-------------------------------------------------------------------
281 :
282 1536 : use tracer_data, only : advance_trcdata
283 : use physics_types,only : physics_state
284 : use ppgrid, only : begchunk, endchunk
285 : use ppgrid, only : pcols, pver
286 : use string_utils, only : to_lower, GLC
287 : use cam_history, only : outfld
288 : use physconst, only : mwdry ! molecular weight dry air ~ kg/kmole
289 : use physconst, only : boltz ! J/K/molecule
290 : ! C.-C. Chen
291 : ! use physconst, only : gravit, rearth
292 : use phys_grid, only : get_wght_all_p
293 :
294 : use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk
295 :
296 : implicit none
297 :
298 : type(physics_state), intent(in) :: state(begchunk:endchunk)
299 :
300 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
301 :
302 370944 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
303 : integer :: ind,c,ncol,i,caseid,m,n
304 : real(r8) :: to_mmr(pcols,pver)
305 370944 : real(r8),pointer :: tmpptr(:,:)
306 :
307 : ! C.-C. Chen
308 : real(r8) :: wght(pcols)
309 :
310 : !------------------------------------------------------------------
311 : ! Return if aircraft_cnt is zero (no aircraft data to process)
312 : !------------------------------------------------------------------
313 370944 : if (aircraft_cnt == 0 ) return
314 0 : call t_startf('All_aircraft_emit_adv')
315 :
316 : !-------------------------------------------------------------------
317 : ! For each field, read more data if needed and interpolate it to the current model time
318 : !-------------------------------------------------------------------
319 0 : do m = 1, aircraft_cnt
320 0 : call advance_trcdata( forcings_air(m)%fields, forcings_air(m)%file, state, pbuf2d)
321 :
322 : !-------------------------------------------------------------------
323 : ! set the tracer fields with the correct units
324 : !-------------------------------------------------------------------
325 0 : do i = 1,number_flds
326 :
327 : ! C.-C. Chen, adding case 4 for kg/sec
328 0 : select case ( to_lower(trim(forcings_air(m)%fields(i)%units(:GLC(forcings_air(m)%fields(i)%units)))) )
329 : case ("molec/cm3","/cm3","molecules/cm3","cm^-3","cm**-3")
330 0 : caseid = 1
331 : case ('kg/kg','mmr')
332 0 : caseid = 2
333 : case ('mol/mol','mole/mole','vmr','fraction')
334 0 : caseid = 3
335 : case ('kg/kg/sec')
336 0 : caseid = 4
337 : case ('kg m-2 s-1')
338 0 : caseid = 5
339 : case ('m/sec' )
340 0 : caseid = 6
341 : case default
342 0 : print*, 'aircraft_emit_adv: units = ',trim(forcings_air(m)%fields(i)%units) ,' are not recognized'
343 0 : call endrun('aircraft_emit_adv: units are not recognized')
344 : end select
345 :
346 0 : ind = index_map(i)
347 :
348 : !$OMP PARALLEL DO PRIVATE (C, NCOL, TO_MMR, tmpptr, pbuf_chnk)
349 0 : do c = begchunk,endchunk
350 0 : ncol = state(c)%ncol
351 :
352 : ! C.-C. Chen, turning emission data to mixing ratio
353 0 : call get_wght_all_p(c,ncol,wght(:ncol))
354 :
355 0 : if (caseid == 1) then
356 0 : to_mmr(:ncol,:) = (molmass(ind)*1.e6_r8*boltz*state(c)%t(:ncol,:))/(mwdry*state(c)%pmiddry(:ncol,:))
357 0 : elseif (caseid == 2) then
358 0 : to_mmr(:ncol,:) = 1._r8
359 0 : elseif (caseid == 4) then
360 : ! do n=1,ncol
361 : ! to_mmr(n,:) = 1.0_r8/(rearth*rearth*wght(n)*state(c)%pdel(n,:)/gravit)
362 : ! end do
363 0 : to_mmr(:ncol,:) = 1.0_r8
364 0 : elseif (caseid == 5) then
365 0 : to_mmr(:ncol,:) = 1.0_r8
366 0 : elseif (caseid == 6) then
367 0 : to_mmr(:ncol,:) = 1.0_r8
368 : else
369 0 : to_mmr(:ncol,:) = molmass(ind)/mwdry
370 : endif
371 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, c)
372 0 : call pbuf_get_field(pbuf_chnk, forcings_air(m)%fields(i)%pbuf_ndx, tmpptr )
373 :
374 0 : tmpptr(:ncol,:) = tmpptr(:ncol,:)*to_mmr(:ncol,:)
375 :
376 0 : call outfld( forcings_air(m)%fields(i)%fldnam, &
377 0 : tmpptr(:ncol,:), ncol, state(c)%lchnk )
378 :
379 : enddo
380 : enddo
381 : enddo
382 :
383 0 : call t_stopf('All_aircraft_emit_adv')
384 370944 : end subroutine aircraft_emit_adv
385 :
386 1536 : subroutine aircraft_emit_readnl(nlfile)
387 : !-------------------------------------------------------------------
388 : ! **** Read in the aircraft_emit namelist *****
389 : !-------------------------------------------------------------------
390 370944 : use namelist_utils, only: find_group_name
391 : use units, only: getunit, freeunit
392 : use mpishorthand
393 :
394 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
395 :
396 : ! Local variables
397 : integer :: unitn, ierr
398 : character(len=*), parameter :: subname = 'aircraft_emit_readnl'
399 :
400 : character(len=256) :: aircraft_specifier(N_AERO)
401 : character(len=256) :: aircraft_datapath
402 : character(len=24) :: aircraft_type
403 :
404 : namelist /aircraft_emit_nl/ aircraft_specifier, aircraft_datapath, aircraft_type
405 : !-----------------------------------------------------------------------------
406 :
407 : ! Initialize namelist variables from local module variables.
408 21504 : aircraft_specifier= air_specifier
409 1536 : aircraft_datapath = air_datapath
410 1536 : aircraft_type = air_type
411 :
412 : ! Read namelist
413 1536 : if (masterproc) then
414 2 : unitn = getunit()
415 2 : open( unitn, file=trim(nlfile), status='old' )
416 2 : call find_group_name(unitn, 'aircraft_emit_nl', status=ierr)
417 2 : if (ierr == 0) then
418 0 : read(unitn, aircraft_emit_nl, iostat=ierr)
419 0 : if (ierr /= 0) then
420 0 : call endrun(subname // ':: ERROR reading namelist')
421 : end if
422 : end if
423 2 : close(unitn)
424 2 : call freeunit(unitn)
425 : end if
426 :
427 : #ifdef SPMD
428 : ! Broadcast namelist variables
429 1536 : call mpibcast(aircraft_specifier,len(aircraft_specifier(1))*N_AERO, mpichar, 0, mpicom)
430 1536 : call mpibcast(aircraft_datapath, len(aircraft_datapath), mpichar, 0, mpicom)
431 1536 : call mpibcast(aircraft_type, len(aircraft_type), mpichar, 0, mpicom)
432 : #endif
433 :
434 : ! Update module variables with user settings.
435 21504 : air_specifier = aircraft_specifier
436 1536 : air_datapath = aircraft_datapath
437 1536 : air_type = aircraft_type
438 :
439 1536 : end subroutine aircraft_emit_readnl
440 :
441 0 : integer function get_aircraft_ndx( name )
442 :
443 : implicit none
444 : character(len=*), intent(in) :: name
445 :
446 : integer :: i
447 :
448 0 : get_aircraft_ndx = 0
449 0 : do i = 1,N_AERO
450 0 : if ( trim(name) == trim(aero_names(i)) ) then
451 0 : get_aircraft_ndx = i
452 0 : return
453 : endif
454 : enddo
455 :
456 : end function get_aircraft_ndx
457 :
458 0 : end module aircraft_emit
|
