Line data Source code
1 : module mo_airplane
2 : !--------------------------------------------------------------------
3 : ! ... Airplane insitu emission sources
4 : !--------------------------------------------------------------------
5 :
6 : use shr_kind_mod, only : r8 => shr_kind_r8
7 : use cam_abortutils, only : endrun
8 : use pio, only : pio_inq_dimid, pio_inq_dimlen, pio_get_var, &
9 : file_desc_t, var_desc_t, pio_inq_vardimid, pio_inq_varndims, pio_nowrite, &
10 : pio_inq_varid, pio_closefile
11 :
12 : use cam_pio_utils, only : cam_pio_openfile
13 : use cam_logfile, only : iulog
14 : implicit none
15 :
16 : private
17 :
18 : save
19 :
20 : real(r8), allocatable :: &
21 : pno(:,:,:), &
22 : pco(:,:,:), &
23 : air_altitude(:)
24 : public :: airpl_set, airpl_src
25 :
26 :
27 :
28 : logical :: has_airpl_src = .false.
29 :
30 : contains
31 :
32 1489176 : subroutine airpl_set( lchnk, ncol, no_ndx, co_ndx, xno_ndx, cldtop, zint_abs, extfrc)
33 : use ppgrid, only : pver
34 : use cam_history, only : outfld
35 :
36 : implicit none
37 : integer, intent(in) :: lchnk, ncol, no_ndx, co_ndx, xno_ndx
38 : real(r8), intent(in) :: cldtop(:), zint_abs(:,:)
39 : real(r8), intent(inout) :: extfrc(:,:,:)
40 :
41 :
42 : ! Local Variables
43 1489176 : real(r8), dimension(ncol,pver) :: no_air, co_air
44 : real(r8) :: ztab_top, ztab_bot, zdel, zdeli, frac, zlev_top, zlev_bot
45 : integer :: nlev
46 : integer :: cldind, kk, i, k
47 :
48 2314006344 : no_air(:,:) = 0._r8
49 2314006344 : co_air(:,:) = 0._r8
50 :
51 1489176 : if(has_airpl_src) then
52 : !---------------------------------------------------------------------
53 : ! ... Add the airplane emissions; must do vertical interpolation
54 : !---------------------------------------------------------------------
55 0 : ztab_top = maxval( air_altitude )
56 0 : ztab_bot = minval( air_altitude )
57 0 : nlev = size(air_altitude) - 1
58 :
59 : !---------------------------------------------------------------------
60 : ! ... add the airplane emissions; must do vertical interpolation
61 : ! Note: the interpolation code is conserving and assumes the
62 : ! aircraft emission vertical grid is uniform with a
63 : ! one kilometer spacing
64 : !---------------------------------------------------------------------
65 0 : level_loop : do k = 1,pver
66 0 : long_loop : do i = 1,ncol
67 0 : zlev_top = zint_abs(i,k) ! altitude at top of model level (km)
68 0 : zlev_bot = zint_abs(i,k+1) ! altitude at bottom of model level (km)
69 0 : zdel = (zlev_top - zlev_bot) * 1.e5_r8 ! model level thickness (cm)
70 0 : zdeli = 1._r8/zdel
71 0 : if( zlev_bot <= ztab_top .and. zlev_top >= ztab_bot ) then
72 0 : do kk = 1,nlev
73 0 : if( zlev_bot <= air_altitude(kk+1) .and. zlev_top >= air_altitude(kk) ) then
74 : frac = (min( zlev_top, air_altitude(kk+1) ) - max( zlev_bot, air_altitude(kk) )) &
75 0 : /(air_altitude(kk+1) - air_altitude(kk)) ! *del_alti(kk)
76 0 : if( no_ndx > 0 ) then
77 0 : extfrc(i,k,no_ndx) = extfrc(i,k,no_ndx) + frac * pno(i,kk,lchnk) * zdeli
78 0 : no_air(i,k) = frac * pno(i,kk,lchnk) * zdeli
79 : end if
80 0 : if( xno_ndx > 0 ) then
81 0 : extfrc(i,k,xno_ndx) = extfrc(i,k,xno_ndx) + frac * pno(i,kk,lchnk) * zdeli
82 : end if
83 0 : if( co_ndx > 0 ) then
84 0 : extfrc(i,k,co_ndx) = extfrc(i,k,co_ndx) + frac * pco(i,kk,lchnk) * zdeli
85 0 : co_air(i,k) = frac * pco(i,kk,lchnk) * zdeli
86 : end if
87 : end if
88 : end do
89 : end if
90 0 : if( k == pver ) then
91 0 : do kk = 1,nlev
92 0 : if( zlev_bot > air_altitude(kk) ) then
93 0 : frac = (min( zlev_bot, air_altitude(kk+1) ) - air_altitude(kk)) &
94 0 : /(air_altitude(kk+1) - air_altitude(kk)) ! *del_alti(kk)
95 0 : if( no_ndx > 0 ) then
96 0 : extfrc(i,k,no_ndx) = extfrc(i,k,no_ndx) + frac * pno(i,kk,lchnk) * zdeli
97 0 : no_air(i,k) = frac * pno(i,kk,lchnk) * zdeli
98 : end if
99 0 : if( xno_ndx > 0 ) then
100 0 : extfrc(i,k,xno_ndx) = extfrc(i,k,xno_ndx) + frac * pno(i,kk,lchnk) * zdeli
101 : end if
102 0 : if( co_ndx > 0 ) then
103 0 : extfrc(i,k,co_ndx) = extfrc(i,k,co_ndx) + frac * pco(i,kk,lchnk) * zdeli
104 0 : co_air(i,k) = frac * pco(i,kk,lchnk) * zdeli
105 : end if
106 : else
107 : exit
108 : end if
109 : end do
110 : end if
111 : end do long_loop
112 : end do level_loop
113 :
114 : end if
115 1489176 : call outfld( 'NO_Aircraft', no_air(:ncol,:), ncol, lchnk )
116 1489176 : call outfld( 'CO_Aircraft', co_air(:ncol,:), ncol, lchnk )
117 :
118 1489176 : end subroutine airpl_set
119 :
120 :
121 1536 : subroutine airpl_src( airpl_emis_file )
122 : !-----------------------------------------------------------------------
123 : ! ... Initialize airplane emissions
124 : ! Note: The emissions are read in in units of molecules/cm**2/s
125 : ! on a vertically resolved grid.
126 : ! Conversion to units of molec/cm**3/s is done in SETEXT
127 : !-----------------------------------------------------------------------
128 1489176 : use spmd_utils, only : masterproc
129 : use interpolate_data, only : lininterp_init, lininterp, lininterp_finish, &
130 : interp_type
131 : use chem_mods, only : adv_mass
132 : use ioFileMod, only : getfil
133 : use mo_chem_utls, only : get_spc_ndx, get_extfrc_ndx
134 : use phys_grid, only : get_ncols_p, get_rlat_all_p, get_rlon_all_p
135 : use ppgrid, only : begchunk, endchunk, pcols
136 : use mo_constants, only : pi, d2r, rearth
137 : use gmean_mod, only : gmean
138 : implicit none
139 :
140 : !-----------------------------------------------------------------------
141 : ! ... Dummy args
142 : !-----------------------------------------------------------------------
143 : character(len=*), intent(in) :: airpl_emis_file
144 :
145 : !-----------------------------------------------------------------------
146 : ! ... Local variables
147 : !-----------------------------------------------------------------------
148 : real(r8), parameter :: msq2cmsq = 1.e4_r8, zero=0._r8, twopi=2._r8*pi
149 : integer :: k
150 : integer :: nlat, nlon, nlev, ndims
151 : integer :: ierr
152 : type(file_desc_t) :: piofile
153 : type(var_desc_t) :: vid
154 : integer :: dimid_lat, dimid_lon, dimid_lev
155 : integer :: dimid(3)
156 1536 : real(r8), allocatable :: lat(:), lon(:)
157 1536 : real(r8), allocatable :: pno_in(:,:,:), pco_in(:,:,:)
158 : real(r8) :: total(2), tmp
159 : real(r8) :: factor
160 : character(len=256) :: locfn
161 : integer :: co_ndx, no_ndx
162 : type(interp_type) :: lon_wgts, lat_wgts
163 : real(r8) :: to_lats(pcols), to_lons(pcols)
164 : integer :: ncols, c
165 :
166 1536 : co_ndx = get_extfrc_ndx('CO')
167 1536 : no_ndx = get_extfrc_ndx('NO')
168 :
169 1536 : if ( co_ndx < 0 .and. no_ndx < 0 ) then
170 1536 : if( masterproc ) then
171 2 : write(iulog,*) 'airpl_src: NO and CO do not have external source --> no aircraft sources will be applied'
172 : endif
173 1536 : return
174 : endif
175 :
176 0 : if ( len_trim(airpl_emis_file) == 0 ) then
177 : return
178 : endif
179 :
180 0 : has_airpl_src = .true.
181 :
182 0 : co_ndx = get_spc_ndx('CO')
183 0 : no_ndx = get_spc_ndx('NO')
184 :
185 :
186 : !-----------------------------------------------------------------------
187 : ! ... Open NetCDF file
188 : !-----------------------------------------------------------------------
189 0 : call getfil (airpl_emis_file, locfn, 0)
190 0 : call cam_pio_openfile (piofile, trim(locfn), PIO_NOWRITE)
191 :
192 : !-----------------------------------------------------------------------
193 : ! ... Get grid dimensions from file
194 : !-----------------------------------------------------------------------
195 0 : ierr = pio_inq_dimid( piofile, 'lat', dimid_lat )
196 0 : ierr = pio_inq_dimlen( piofile, dimid_lat, nlat )
197 0 : allocate( lat(nlat), stat=ierr )
198 0 : if( ierr /= 0 ) then
199 0 : write(iulog,*) 'airpl_src: lat allocation error = ',ierr
200 0 : call endrun
201 : end if
202 0 : ierr = pio_inq_varid( piofile, 'lat', vid )
203 0 : ierr = pio_get_var( piofile, vid, lat )
204 0 : lat(:nlat) = lat(:nlat) * d2r
205 :
206 0 : ierr = pio_inq_dimid( piofile, 'lon', dimid_lon )
207 0 : ierr = pio_inq_dimlen( piofile, dimid_lon, nlon )
208 0 : allocate( lon(nlon), stat=ierr )
209 0 : if( ierr /= 0 ) then
210 0 : write(iulog,*) 'airpl_src: lon allocation error = ',ierr
211 0 : call endrun
212 : end if
213 0 : ierr = pio_inq_varid( piofile, 'lon', vid )
214 0 : ierr = pio_get_var( piofile, vid, lon )
215 0 : lon(:nlon) = lon(:nlon) * d2r
216 :
217 0 : ierr = pio_inq_dimid( piofile, 'altitude', dimid_lev )
218 0 : ierr = pio_inq_dimlen( piofile, dimid_lev, nlev )
219 0 : allocate( air_altitude(nlev+1), stat=ierr )
220 0 : if( ierr /= 0 ) then
221 0 : write(iulog,*) 'airpl_src: air_altitude allocation error = ',ierr
222 0 : call endrun
223 : end if
224 0 : ierr = pio_inq_varid( piofile, 'altitude', vid )
225 0 : ierr = pio_get_var( piofile, vid, air_altitude(1:nlev) )
226 0 : air_altitude(nlev+1) = air_altitude(nlev) + (air_altitude(nlev) - air_altitude(nlev-1))
227 :
228 : !-----------------------------------------------------------------------
229 : ! ... Set up regridding
230 : !-----------------------------------------------------------------------
231 :
232 0 : allocate( pno_in(nlon,nlat,nlev), stat=ierr )
233 0 : if( ierr /= 0 ) then
234 0 : write(iulog,*) 'airpl_src: pno_in allocation error = ',ierr
235 0 : call endrun
236 : end if
237 0 : allocate( pco_in(nlon,nlat,nlev), stat=ierr )
238 0 : if( ierr /= 0 ) then
239 0 : write(iulog,*) 'airpl_src: pco_in allocation error = ',ierr
240 0 : call endrun
241 : end if
242 0 : allocate(pno(pcols,nlev,begchunk:endchunk), stat=ierr )
243 0 : if( ierr /= 0 ) then
244 0 : write(iulog,*) 'airpl_src: pno allocation error = ',ierr
245 0 : call endrun
246 : end if
247 0 : allocate( pco(pcols,nlev,begchunk:endchunk), stat=ierr )
248 0 : if( ierr /= 0 ) then
249 0 : write(iulog,*) 'airpl_src: pco allocation error = ',ierr
250 0 : call endrun
251 : end if
252 :
253 : !-----------------------------------------------------------------------
254 : ! ... Read emissions
255 : !-----------------------------------------------------------------------
256 0 : ierr = pio_inq_varid( piofile, 'nox', vid )
257 0 : ierr = pio_inq_varndims( piofile, vid, ndims )
258 0 : if( ndims /= 3 ) then
259 0 : write(iulog,*) 'airpl_src: variable nox has ndims = ',ndims,', expecting 3'
260 0 : call endrun
261 : end if
262 0 : ierr = pio_inq_vardimid( piofile, vid, dimid )
263 0 : if( dimid(1) /= dimid_lon .or. dimid(2) /= dimid_lat .or. dimid(3) /= dimid_lev ) then
264 0 : write(iulog,*) 'airpl_src: Dimensions in wrong order for variable nox'
265 0 : write(iulog,*) '... Expecting (lon, lat, lev)'
266 0 : call endrun
267 : end if
268 : ierr = pio_get_var( piofile, vid, &
269 : (/ 1, 1, 1/), & ! start
270 : (/ nlon, nlat, nlev /), & ! count
271 0 : pno_in )
272 :
273 0 : ierr = pio_inq_varid( piofile, 'co', vid )
274 0 : ierr = pio_inq_varndims( piofile, vid, ndims )
275 :
276 0 : if( ndims /= 3 ) then
277 0 : write(iulog,*) 'READ_SFLX: variable co has ndims = ',ndims,', expecting 3'
278 0 : call endrun
279 : end if
280 0 : ierr = pio_inq_vardimid( piofile, vid, dimid )
281 0 : if( dimid(1) /= dimid_lon .or. dimid(2) /= dimid_lat .or. dimid(3) /= dimid_lev ) then
282 0 : write(iulog,*) 'airpl_src: Dimensions in wrong order for variable co'
283 0 : write(iulog,*) '... Expecting (lon, lat, lev)'
284 0 : call endrun
285 : end if
286 : ierr = pio_get_var( piofile, vid, &
287 : (/ 1, 1, 1/), & ! start
288 : (/ nlon, nlat, nlev /), & ! count
289 0 : pco_in )
290 0 : call pio_closefile( piofile )
291 :
292 : !-----------------------------------------------------------------------
293 : ! ... Regrid emissions
294 : !-----------------------------------------------------------------------
295 0 : do c=begchunk,endchunk
296 0 : ncols = get_ncols_p(c)
297 0 : call get_rlat_all_p(c, pcols, to_lats)
298 0 : call get_rlon_all_p(c, pcols, to_lons)
299 0 : call lininterp_init(lon, nlon, to_lons, ncols, 2, lon_wgts, zero, twopi)
300 0 : call lininterp_init(lat, nlat, to_lats, ncols, 1, lat_wgts)
301 :
302 0 : do k = 1,nlev
303 0 : call lininterp(pno_in(:,:,k), nlon, nlat, pno(:,k,c), ncols, lon_wgts, lat_wgts)
304 0 : call lininterp(pco_in(:,:,k), nlon, nlat, pco(:,k,c), ncols, lon_wgts, lat_wgts)
305 : enddo
306 0 : call lininterp_finish(lon_wgts)
307 0 : call lininterp_finish(lat_wgts)
308 : enddo
309 :
310 0 : deallocate( pno_in, pco_in, lon, lat, stat=ierr )
311 0 : if( ierr /= 0 ) then
312 0 : write(iulog,*) 'airpl_src: Failed to deallocate pno_in,pco_in; ierr = ',ierr
313 0 : call endrun
314 : end if
315 : !-----------------------------------------------------------------------
316 : ! ... Get global emission from this source
317 : !-----------------------------------------------------------------------
318 0 : total = zero
319 0 : do k=1,nlev
320 0 : call gmean(pno(:,k,:), tmp)
321 0 : total(1)= total(1)+tmp
322 0 : call gmean(pco(:,k,:), tmp)
323 0 : total(2)= total(2)+tmp
324 : end do
325 :
326 0 : if(masterproc) then
327 :
328 : factor = 86400._r8 * 365._r8 & ! sec / year
329 : / 6.022e23_r8 & ! molec / mole
330 : * 1.e-12_r8 & ! Tg / g
331 : * msq2cmsq & ! meters**2 to cm**2
332 0 : * 4._r8*pi*rearth*rearth ! global mean to global total
333 :
334 0 : write(iulog,*) 'airpl_src: nlev = ',nlev
335 : !-----------------------------------------------------------------------
336 : ! ... Convert totals from molec cm^-2 s^-1 to Tg y^-1
337 : !-----------------------------------------------------------------------
338 0 : if (no_ndx .gt. 0) then
339 0 : total(1) = total(1) * adv_mass(no_ndx) * factor
340 0 : write(iulog,'('' airpl_src Aircraft emissions: '',a6,'' = '',f10.3,1X,a6)') 'NO',total(1),'TgN/y'
341 : endif
342 0 : if (co_ndx .gt. 0) then
343 0 : total(2) = total(2) * adv_mass(co_ndx) * factor
344 0 : write(iulog,'('' airpl_src Aircraft emissions: '',a6,'' = '',f10.3,1X,a6)') 'CO',total(2),'Tg/y'
345 : endif
346 : end if
347 :
348 3072 : end subroutine airpl_src
349 :
350 : end module mo_airplane
|
