Line data Source code
1 : module mo_srf_emissions
2 : !---------------------------------------------------------------
3 : ! ... surface emissions module
4 : !---------------------------------------------------------------
5 :
6 : use shr_kind_mod, only : r8 => shr_kind_r8
7 : use chem_mods, only : gas_pcnst
8 : use spmd_utils, only : masterproc
9 : use cam_abortutils,only : endrun
10 : use ioFileMod, only : getfil
11 : use cam_logfile, only : iulog
12 : use tracer_data, only : trfld,trfile
13 :
14 : implicit none
15 :
16 : type :: emission
17 : integer :: spc_ndx
18 : real(r8) :: mw
19 : real(r8) :: scalefactor
20 : character(len=256):: filename
21 : character(len=16) :: species
22 : character(len=8) :: units
23 : integer :: nsectors
24 : character(len=32),pointer :: sectors(:)
25 : type(trfld), pointer :: fields(:)
26 : type(trfile) :: file
27 : end type emission
28 :
29 : private
30 :
31 : public :: srf_emissions_inti, set_srf_emissions, set_srf_emissions_time
32 :
33 : logical, public, protected :: has_emis(gas_pcnst) = .false.
34 :
35 : real(r8), parameter :: amufac = 1.65979e-23_r8 ! 1.e4* kg / amu
36 : type(emission), allocatable :: emissions(:)
37 : integer :: n_emis_files
38 : integer :: c10h16_ndx, isop_ndx
39 :
40 : contains
41 :
42 1536 : subroutine srf_emissions_inti( srf_emis_specifier, emis_type_in, emis_cycle_yr, emis_fixed_ymd, emis_fixed_tod )
43 :
44 : !-----------------------------------------------------------------------
45 : ! ... initialize the surface emissions
46 : !-----------------------------------------------------------------------
47 :
48 : use chem_mods, only : adv_mass
49 : use mo_chem_utls, only : get_spc_ndx
50 : use tracer_data, only : trcdata_init
51 : use cam_pio_utils, only : cam_pio_openfile
52 : use pio, only : pio_inquire, pio_nowrite, pio_closefile, pio_inq_varndims
53 : use pio, only : pio_inq_varname, pio_inq_vardimid, pio_inq_dimid
54 : use pio, only : file_desc_t, pio_get_att, PIO_NOERR, PIO_GLOBAL
55 : use pio, only : pio_seterrorhandling, PIO_BCAST_ERROR,PIO_INTERNAL_ERROR
56 : use chem_surfvals, only : flbc_list
57 : use string_utils, only : GLC
58 : use m_MergeSorts, only : IndexSort
59 :
60 : implicit none
61 :
62 : !-----------------------------------------------------------------------
63 : ! ... dummy arguments
64 : !-----------------------------------------------------------------------
65 : character(len=*), intent(in) :: srf_emis_specifier(:)
66 : character(len=*), intent(in) :: emis_type_in
67 : integer, intent(in) :: emis_cycle_yr
68 : integer, intent(in) :: emis_fixed_ymd
69 : integer, intent(in) :: emis_fixed_tod
70 :
71 : !-----------------------------------------------------------------------
72 : ! ... local variables
73 : !-----------------------------------------------------------------------
74 : integer :: astat
75 : integer :: j, l, m, n, i, nn ! Indices
76 : character(len=16) :: spc_name
77 : character(len=256) :: filename
78 :
79 3072 : character(len=16) :: emis_species(size(srf_emis_specifier))
80 3072 : character(len=256) :: emis_filenam(size(srf_emis_specifier))
81 3072 : integer :: emis_indexes(size(srf_emis_specifier))
82 3072 : integer :: indx(size(srf_emis_specifier))
83 3072 : real(r8) :: emis_scalefactor(size(srf_emis_specifier))
84 :
85 : integer :: vid, nvars, isec, num_dims_emis
86 : integer :: vndims
87 1536 : logical, allocatable :: is_sector(:)
88 : type(file_desc_t) :: ncid
89 : character(len=32) :: varname
90 : character(len=256) :: locfn
91 : integer :: ierr
92 : character(len=1), parameter :: filelist = ''
93 : character(len=1), parameter :: datapath = ''
94 : logical , parameter :: rmv_file = .false.
95 : logical :: unstructured
96 : character(len=32) :: emis_type = ' '
97 : character(len=80) :: file_interp_type = ' '
98 : character(len=256) :: tmp_string = ' '
99 : character(len=32) :: xchr = ' '
100 : real(r8) :: xdbl
101 : integer :: time_dimid, ncol_dimid
102 1536 : integer, allocatable :: dimids(:)
103 :
104 1536 : has_emis(:) = .false.
105 1536 : nn = 0
106 155136 : indx(:) = 0
107 :
108 47616 : count_emis: do n=1,size(srf_emis_specifier)
109 47616 : if ( len_trim(srf_emis_specifier(n) ) == 0 ) then
110 : exit count_emis
111 : endif
112 :
113 46080 : i = scan(srf_emis_specifier(n),'->')
114 46080 : spc_name = trim(adjustl(srf_emis_specifier(n)(:i-1)))
115 :
116 : ! need to parse out scalefactor ...
117 46080 : tmp_string = adjustl(srf_emis_specifier(n)(i+2:))
118 46080 : j = scan( tmp_string, '*' )
119 46080 : if (j>0) then
120 18432 : xchr = tmp_string(1:j-1) ! get the multipler (left of the '*')
121 18432 : read( xchr, * ) xdbl ! convert the string to a real
122 18432 : tmp_string = adjustl(tmp_string(j+1:)) ! get the filepath name (right of the '*')
123 : else
124 27648 : xdbl = 1._r8
125 : endif
126 46080 : filename = trim(tmp_string)
127 :
128 46080 : m = get_spc_ndx(spc_name)
129 :
130 46080 : if (m > 0) then
131 46080 : has_emis(m) = .true.
132 : else
133 0 : write(iulog,*) 'srf_emis_inti: spc_name ',spc_name,' is not included in the simulation'
134 0 : call endrun('srf_emis_inti: invalid surface emission specification')
135 : endif
136 :
137 1935360 : if (any( flbc_list == spc_name )) then
138 : call endrun('srf_emis_inti: ERROR -- cannot specify both fixed LBC ' &
139 0 : //'and emissions for the same species: '//trim(spc_name))
140 : endif
141 :
142 46080 : nn = nn+1
143 46080 : emis_species(nn) = spc_name
144 46080 : emis_filenam(nn) = filename
145 46080 : emis_indexes(nn) = m
146 46080 : emis_scalefactor(nn) = xdbl
147 :
148 47616 : indx(n)=n
149 :
150 : enddo count_emis
151 :
152 1536 : n_emis_files = nn
153 :
154 1536 : if (masterproc) write(iulog,*) 'srf_emis_inti: n_emis_files = ',n_emis_files
155 :
156 56832 : allocate( emissions(n_emis_files), stat=astat )
157 1536 : if( astat/= 0 ) then
158 0 : write(iulog,*) 'srf_emis_inti: failed to allocate emissions array; error = ',astat
159 0 : call endrun('srf_emis_inti: failed to allocate emissions array')
160 : end if
161 :
162 : !-----------------------------------------------------------------------
163 : ! Sort the input files so that the emissions sources are summed in the
164 : ! same order regardless of the order of the input files in the namelist
165 : !-----------------------------------------------------------------------
166 1536 : if (n_emis_files > 0) then
167 1536 : call IndexSort(n_emis_files, indx, emis_filenam)
168 : end if
169 :
170 : !-----------------------------------------------------------------------
171 : ! ... setup the emission type array
172 : !-----------------------------------------------------------------------
173 47616 : do m=1,n_emis_files
174 46080 : emissions(m)%spc_ndx = emis_indexes(indx(m))
175 46080 : emissions(m)%units = 'Tg/y'
176 46080 : emissions(m)%species = emis_species(indx(m))
177 46080 : emissions(m)%mw = adv_mass(emis_indexes(indx(m))) ! g / mole
178 46080 : emissions(m)%filename = emis_filenam(indx(m))
179 47616 : emissions(m)%scalefactor = emis_scalefactor(indx(m))
180 : enddo
181 :
182 : !-----------------------------------------------------------------------
183 : ! read emis files to determine number of sectors
184 : !-----------------------------------------------------------------------
185 47616 : spc_loop: do m = 1, n_emis_files
186 :
187 46080 : emissions(m)%nsectors = 0
188 :
189 46080 : if (masterproc) then
190 60 : write(iulog,'(a,i3,a)') 'srf_emissions_inti m: ',m,' init file : '//trim(emissions(m)%filename)
191 : endif
192 :
193 46080 : call getfil (emissions(m)%filename, locfn, 0)
194 46080 : call cam_pio_openfile ( ncid, trim(locfn), PIO_NOWRITE)
195 46080 : ierr = pio_inquire (ncid, nVariables=nvars)
196 :
197 46080 : call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
198 46080 : ierr = pio_inq_dimid( ncid, 'ncol', ncol_dimid )
199 46080 : unstructured = ierr==PIO_NOERR
200 46080 : call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
201 :
202 138240 : allocate(is_sector(nvars))
203 328704 : is_sector(:) = .false.
204 :
205 46080 : if (unstructured) then
206 46080 : ierr = pio_inq_dimid( ncid, 'time', time_dimid )
207 : end if
208 :
209 328704 : do vid = 1,nvars
210 :
211 282624 : ierr = pio_inq_varndims (ncid, vid, vndims)
212 :
213 282624 : if (unstructured) then
214 : num_dims_emis = 2
215 : else
216 0 : num_dims_emis = 3
217 : endif
218 :
219 282624 : if( vndims < num_dims_emis ) then
220 : cycle
221 52224 : elseif( vndims > num_dims_emis ) then
222 0 : ierr = pio_inq_varname (ncid, vid, varname)
223 0 : write(iulog,*) 'srf_emis_inti: Skipping variable ', trim(varname),', ndims = ',vndims, &
224 0 : ' , species=',trim(emissions(m)%species)
225 0 : cycle
226 : end if
227 :
228 98304 : if (unstructured) then
229 156672 : allocate( dimids(vndims) )
230 52224 : ierr = pio_inq_vardimid( ncid, vid, dimids )
231 104448 : if ( any(dimids(:)==ncol_dimid) .and. any(dimids(:)==time_dimid) ) then
232 52224 : emissions(m)%nsectors = emissions(m)%nsectors+1
233 52224 : is_sector(vid)=.true.
234 : endif
235 52224 : deallocate(dimids)
236 : else
237 0 : emissions(m)%nsectors = emissions(m)%nsectors+1
238 0 : is_sector(vid)=.true.
239 : end if
240 :
241 : enddo
242 :
243 138240 : allocate( emissions(m)%sectors(emissions(m)%nsectors), stat=astat )
244 46080 : if( astat/= 0 ) then
245 0 : write(iulog,*) 'srf_emis_inti: failed to allocate emissions(m)%sectors array; error = ',astat
246 0 : call endrun
247 : end if
248 :
249 46080 : isec = 1
250 :
251 328704 : do vid = 1,nvars
252 328704 : if( is_sector(vid) ) then
253 52224 : ierr = pio_inq_varname(ncid, vid, emissions(m)%sectors(isec))
254 52224 : isec = isec+1
255 : endif
256 : enddo
257 46080 : deallocate(is_sector)
258 :
259 : ! Global attribute 'input_method' overrides the srf_emis_type namelist setting on
260 : ! a file-by-file basis. If the emis file does not contain the 'input_method'
261 : ! attribute then the srf_emis_type namelist setting is used.
262 46080 : call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
263 46080 : ierr = pio_get_att(ncid, PIO_GLOBAL, 'input_method', file_interp_type)
264 46080 : call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
265 46080 : if ( ierr == PIO_NOERR) then
266 0 : l = GLC(file_interp_type)
267 0 : emis_type(1:l) = file_interp_type(1:l)
268 0 : emis_type(l+1:) = ' '
269 : else
270 46080 : emis_type = trim(emis_type_in)
271 : endif
272 :
273 46080 : call pio_closefile (ncid)
274 :
275 138240 : allocate(emissions(m)%file%in_pbuf(size(emissions(m)%sectors)))
276 98304 : emissions(m)%file%in_pbuf(:) = .false.
277 :
278 0 : call trcdata_init( emissions(m)%sectors, &
279 : emissions(m)%filename, filelist, datapath, &
280 : emissions(m)%fields, &
281 : emissions(m)%file, &
282 47616 : rmv_file, emis_cycle_yr, emis_fixed_ymd, emis_fixed_tod, trim(emis_type) )
283 :
284 : enddo spc_loop
285 :
286 1536 : c10h16_ndx = get_spc_ndx('C10H16')
287 1536 : isop_ndx = get_spc_ndx('ISOP')
288 :
289 3072 : end subroutine srf_emissions_inti
290 :
291 741888 : subroutine set_srf_emissions_time( pbuf2d, state )
292 : !-----------------------------------------------------------------------
293 : ! ... check serial case for time span
294 : !-----------------------------------------------------------------------
295 :
296 1536 : use physics_types,only : physics_state
297 : use ppgrid, only : begchunk, endchunk
298 : use tracer_data, only : advance_trcdata
299 : use physics_buffer, only : physics_buffer_desc
300 :
301 : implicit none
302 :
303 : type(physics_state), intent(in):: state(begchunk:endchunk)
304 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
305 :
306 : !-----------------------------------------------------------------------
307 : ! ... local variables
308 : !-----------------------------------------------------------------------
309 : integer :: m
310 :
311 11499264 : do m = 1,n_emis_files
312 11499264 : call advance_trcdata( emissions(m)%fields, emissions(m)%file, state, pbuf2d )
313 : end do
314 :
315 370944 : end subroutine set_srf_emissions_time
316 :
317 : ! adds surf flux specified in file to sflx
318 1489176 : subroutine set_srf_emissions( lchnk, ncol, sflx )
319 : !--------------------------------------------------------
320 : ! ... form the surface fluxes for this latitude slice
321 : !--------------------------------------------------------
322 :
323 370944 : use mo_constants, only : pi
324 : use time_manager, only : get_curr_calday
325 : use string_utils, only : to_lower, GLC
326 : use phys_grid, only : get_rlat_all_p, get_rlon_all_p
327 :
328 : implicit none
329 :
330 : !--------------------------------------------------------
331 : ! ... Dummy arguments
332 : !--------------------------------------------------------
333 : integer, intent(in) :: ncol ! columns in chunk
334 : integer, intent(in) :: lchnk ! chunk index
335 : real(r8), intent(out) :: sflx(:,:) ! surface emissions ( kg/m^2/s )
336 :
337 : !--------------------------------------------------------
338 : ! ... local variables
339 : !--------------------------------------------------------
340 : integer :: i, m, n
341 : real(r8) :: factor
342 : real(r8) :: dayfrac ! fration of day in light
343 : real(r8) :: iso_off ! time iso flux turns off
344 : real(r8) :: iso_on ! time iso flux turns on
345 :
346 : logical :: polar_day,polar_night
347 : real(r8) :: doy_loc
348 : real(r8) :: sunon,sunoff
349 : real(r8) :: loc_angle
350 : real(r8) :: latitude
351 : real(r8) :: declination
352 : real(r8) :: tod
353 : real(r8) :: calday
354 :
355 : real(r8), parameter :: dayspy = 365._r8
356 : real(r8), parameter :: twopi = 2.0_r8 * pi
357 : real(r8), parameter :: pid2 = 0.5_r8 * pi
358 : real(r8), parameter :: dec_max = 23.45_r8 * pi/180._r8
359 :
360 2978352 : real(r8) :: flux(ncol)
361 : real(r8) :: mfactor
362 : integer :: isec
363 :
364 : character(len=12),parameter :: mks_units(4) = (/ "kg/m2/s ", &
365 : "kg/m2/sec ", &
366 : "kg/m^2/s ", &
367 : "kg/m^2/sec " /)
368 : character(len=12) :: units
369 :
370 2978352 : real(r8), dimension(ncol) :: rlats, rlons
371 :
372 786284928 : sflx(:,:) = 0._r8
373 :
374 : !--------------------------------------------------------
375 : ! ... set non-zero emissions
376 : !--------------------------------------------------------
377 46164456 : emis_loop : do m = 1,n_emis_files
378 :
379 44675280 : n = emissions(m)%spc_ndx
380 :
381 745973280 : flux(:) = 0._r8
382 95307264 : do isec = 1,emissions(m)%nsectors
383 890111664 : flux(:ncol) = flux(:ncol) + emissions(m)%scalefactor*emissions(m)%fields(isec)%data(:ncol,1,lchnk)
384 : enddo
385 :
386 44675280 : units = to_lower(trim(emissions(m)%fields(1)%units(:GLC(emissions(m)%fields(1)%units))))
387 :
388 224865576 : if ( any( mks_units(:) == units ) ) then
389 0 : sflx(:ncol,n) = sflx(:ncol,n) + flux(:ncol)
390 : else
391 44675280 : mfactor = amufac * emissions(m)%mw
392 745973280 : sflx(:ncol,n) = sflx(:ncol,n) + flux(:ncol) * mfactor
393 : endif
394 :
395 : end do emis_loop
396 :
397 1489176 : call get_rlat_all_p( lchnk, ncol, rlats )
398 1489176 : call get_rlon_all_p( lchnk, ncol, rlons )
399 :
400 1489176 : calday = get_curr_calday()
401 1489176 : doy_loc = aint( calday )
402 1489176 : declination = dec_max * cos((doy_loc - 172._r8)*twopi/dayspy)
403 1489176 : tod = (calday - doy_loc) + .5_r8
404 :
405 24865776 : do i = 1,ncol
406 : !
407 23376600 : polar_day = .false.
408 23376600 : polar_night = .false.
409 : !
410 23376600 : loc_angle = tod * twopi + rlons(i)
411 23376600 : loc_angle = mod( loc_angle,twopi )
412 23376600 : latitude = rlats(i)
413 : !
414 : !------------------------------------------------------------------
415 : ! determine if in polar day or night
416 : ! if not in polar day or night then
417 : ! calculate terminator longitudes
418 : !------------------------------------------------------------------
419 23376600 : if( abs(latitude) >= (pid2 - abs(declination)) ) then
420 1575552 : if( sign(1._r8,declination) == sign(1._r8,latitude) ) then
421 : polar_day = .true.
422 : sunoff = 2._r8*twopi
423 : sunon = -twopi
424 : else
425 787776 : polar_night = .true.
426 : end if
427 : else
428 21801048 : sunoff = acos( -tan(declination)*tan(latitude) )
429 21801048 : sunon = twopi - sunoff
430 : end if
431 :
432 : !--------------------------------------------------------
433 : ! ... adjust alpha-pinene for diurnal variation
434 : !--------------------------------------------------------
435 23376600 : if( c10h16_ndx > 0 ) then
436 0 : if( has_emis(c10h16_ndx) ) then
437 0 : if( .not. polar_night .and. .not. polar_day ) then
438 0 : dayfrac = sunoff / pi
439 0 : sflx(i,c10h16_ndx) = sflx(i,c10h16_ndx) / (.7_r8 + .3_r8*dayfrac)
440 0 : if( loc_angle >= sunoff .and. loc_angle <= sunon ) then
441 0 : sflx(i,c10h16_ndx) = sflx(i,c10h16_ndx) * .7_r8
442 : endif
443 : end if
444 : end if
445 : end if
446 :
447 : !--------------------------------------------------------
448 : ! ... adjust isoprene for diurnal variation
449 : !--------------------------------------------------------
450 24865776 : if( isop_ndx > 0 ) then
451 0 : if( has_emis(isop_ndx) ) then
452 0 : if( .not. polar_night ) then
453 0 : if( polar_day ) then
454 : iso_off = .8_r8 * pi
455 : iso_on = 1.2_r8 * pi
456 : else
457 0 : iso_off = .8_r8 * sunoff
458 0 : iso_on = 2._r8 * pi - iso_off
459 : end if
460 0 : if( loc_angle >= iso_off .and. loc_angle <= iso_on ) then
461 0 : sflx(i,isop_ndx) = 0._r8
462 : else
463 0 : factor = loc_angle - iso_on
464 0 : if( factor <= 0._r8 ) then
465 0 : factor = factor + 2._r8*pi
466 : end if
467 0 : factor = factor / (2._r8*iso_off + 1.e-6_r8)
468 0 : sflx(i,isop_ndx) = sflx(i,isop_ndx) * 2._r8 / iso_off * pi * (sin(pi*factor))**2
469 : end if
470 : else
471 0 : sflx(i,isop_ndx) = 0._r8
472 : end if
473 : end if
474 : end if
475 :
476 : end do
477 :
478 1489176 : end subroutine set_srf_emissions
479 :
480 1489176 : end module mo_srf_emissions
|
