./PaxHeaders.31608/01.prepemis.pro0000644000000000000000000000012713235103020013527 xustar0030 mtime=1517585936.042618664 27 atime=1517585936.041825 30 ctime=1517585936.042618664 01.prepemis.pro0000755001151300017500000002436713235103020014127 0ustar00barthmncar00000000000000; code to put aircraft emissions (1x1deg) into WRF-Chem ; December 2008 (Jeff Lee) ;------------------------------------------------------------------------ FUNCTION read_data, file, month, lonmin, lonmax, latmin, latmax id=ncdf_open(file) ncdf_varget, id, 'alt',alt ncdf_varget, id, 'lon',lon ncdf_varget, id, 'lat',lat ncdf_varget, id, 'NO', e_no ncdf_varget, id, 'CO', e_co ncdf_varget, id, 'CH4',e_ch4 ncdf_varget, id, 'SO2',e_so2 ncdf_close,id ; select only one month imonth = month-1 Print, 'Selected month index (0-11): ',imonth e_no_month = reform (e_no[*,*,*,imonth]) e_co_month = reform (e_co[*,*,*,imonth]) e_ch4_month = reform (e_ch4[*,*,*,imonth]) e_so2_month = reform (e_so2[*,*,*,imonth]) ;help, e_no_month ; change longitude from 0~360 to -180~180 k=where(lon ge 180) if k[0] ne -1 then lon[k] = lon[k]-360 ; reduce from global to region of interest ;klon = where(lon ge 0) ;klat = where(lat ge -90, nklat) ;klon = where(lon le 40 or lon ge 300) klon = where(lon ge lonmin and lon le lonmax, nklon) klat = where(lat ge latmin and lat le latmax, nklat) lon=lon[klon] lat=lat[klat] Print, 'nklon= ', nklon, lon Print, 'nklat= ', nklat, lat e_no_month = e_no_month [klon,klat,*] e_co_month = e_co_month [klon,klat,*] e_ch4_month = e_ch4_month[klon,klat,*] e_so2_month = e_so2_month[klon,klat,*] str={alt:alt,lon:lon,lat:lat,e_no:e_no_month,e_co:e_co_month,e_ch4:e_ch4_month,e_so2:e_so2_month} return, str end ;------------------------------------------------------------------------ PRO mean_height, emis_height_max, ph, phb, hgt, eta_dim, wrfheight, wrflevs, dz ; all units in meter print,emis_height_max alt = (ph+phb)/9.81 for i=0,eta_dim-1 do alt[*,*,i]=alt[*,*,i]-hgt mean_alt = fltarr(eta_dim) for i=0,eta_dim-1 do begin mean_alt[i] = mean(alt[*,*,i]-alt[*,*,0]) ;print,i, mean_alt[i], stddev(alt[*,*,i]-alt[*,*,0]) endfor index = where(mean_alt lt emis_height_max) wrflevs= n_elements(index) print, wrflevs, index wrfheight =fltarr(wrflevs) wrfheight = mean_alt[index] print,'' print,'emis_level model mean_height' for i=0,wrflevs-1 do print,i,wrfheight[i] dz = fltarr(wrflevs) print,'' print,'emis_level dz' for i=0,wrflevs-2 do begin dz[i] = mean_alt[i+1]-mean_alt[i] print,i, dz[i] endfor end ;------------------------------------------------------------------------- PRO distance, lat1,lon1,lat2,lon2, d r = 6.370e6 cutpoint = 1.e-9 degToRad = 0.017453293 rlat1 = lat1 * degToRad rlon1 = -1*lon1 * degToRAD rlat2 = lat2 * degToRad rlon2 =-1*lon2 * degToRad dellat = (rlat2 - rlat1) dellon = (rlon2 - rlon1) if ( abs(dellat) lE cutpoint ) then d=0 if ( abs(dellat) gt cutpoint ) then begin a = sin (dellat/2.) * sin(dellat/2.) + cos(rlat1)*cos(rlat2) $ * sin(dellon/2.) * sin(dellon/2.) if (a lt 1.00) then begin c = 2.*asin(sqrt(a) ) endif else if a ge 1.0 then begin c = 2. * asin(1.00) endif d = r * c endif end @mapcf.pro @ijll_lc.pro @llij_lc.pro @create_netcdffile.pro ; Begin of Main Program ;############################################################################ ;user modified: ;input wrfinput file wrfinput_file = './wrfinput_d01' ;input emission file data_file = './baughcum_1999.aircraft.1x1.nc' minhour=0 ; which hours of the day emissions should be created ; e.g. minhour=0 and maxhour=23 creates emissions for entire day maxhour=23 ;pick the month of emissions ;mm = '06' domain = 'd01' yyyy = '2012' mm = '05' dd = '21' ;set to zero below maxheight, which must be the same as NEI emission maxheight = 1000. ;end of user modified ;############################################################################## output_dir = './' ; get infomation from wrfinput file id = ncdf_open(wrfinput_file) ncdf_attget,id,'DX',dx_temp,/GLOBAL ncdf_attget,id,'WEST-EAST_GRID_DIMENSION',wrf_ii,/GLOBAL ncdf_attget,id,'SOUTH-NORTH_GRID_DIMENSION',wrf_jj,/GLOBAL ncdf_attget,id,'BOTTOM-TOP_GRID_DIMENSION',eta_dim,/GLOBAL ;;; GGP ncdf_attget,id,'CEN_LAT',cen_lat,/GLOBAL ncdf_attget,id,'CEN_LON',cen_lon,/GLOBAL ncdf_attget,id,'STAND_LON',stand_lon,/GLOBAL ncdf_attget,id,'TRUELAT1',truelat1,/GLOBAL ncdf_attget,id,'TRUELAT2',truelat2,/GLOBAL ncdf_attget,id,'MAP_PROJ',map_proj,/GLOBAL ncdf_attget,id,'MOAD_CEN_LAT',moad_cen_lat,/GLOBAL ncdf_varget,id,'PH',ph ncdf_varget,id,'PHB',phb ncdf_varget,id,'HGT',hgt ncdf_varget,id,'XLONG',xlong ncdf_varget,id,'XLAT',xlat ncdf_close,id ; set local values wrf_ix = wrf_ii - 1 wrf_jx = wrf_jj - 1 dx = dx_temp/1000. hemi = 1 lat1 = cen_lat lon1 = cen_lon knowni = (wrf_ix+1)/2. knownj = (wrf_jx+1)/2. ; calculates lat/lon for WRF domain wrf_lat = fltarr(wrf_ix, wrf_jx) wrf_lon = fltarr(wrf_ix, wrf_jx) wrf_lat_ll = fltarr(wrf_ix, wrf_jx) wrf_lon_ll = fltarr(wrf_ix, wrf_jx) for i=0,wrf_ix-1 do begin for j=0,wrf_jx-1 do begin ijll_lc, i+1.5,j+1.5, a, b, truelat1, truelat2, $ hemi, stand_lon,lat1, lon1, knowni, knownj, dx wrf_lat[i,j] = a wrf_lon[i,j] = b ijll_lc, i+1.,j+1., a, b, truelat1, truelat2, $ hemi, stand_lon,lat1, lon1, knowni, knownj, dx wrf_lat_ll[i,j] = a wrf_lon_ll[i,j] = b endfor endfor ; get latmin,latmax,lonmin and lonmax ; rough domain borders (to limit the amount of aircraft data to deal with) ; (pick 5 degree extension on all 4 sides) degree_extention = 5.0 i=0 j=0 ijll_lc, i,j, a, b, truelat1, truelat2, $ hemi, stand_lon,lat1, lon1, knowni, knownj, dx latmin = round(a - degree_extention) lonmin = round(b - degree_extention) ;print, latmin, lonmin i=wrf_ix - 1 j=wrf_jx - 1 ijll_lc, i,j, a, b, truelat1, truelat2, $ hemi, stand_lon,lat1, lon1, knowni, knownj, dx latmax = round(a + degree_extention) lonmax = round(b + degree_extention) ;print, latmax, lonmax ; read emission data data = read_data(data_file, fix(mm), lonmin, lonmax, latmin, latmax ) ; find the top altitude (meter) from emission data ; data.alt in km, change it into meter data_alt = data.alt * 1.e3 print, data_alt n_alt = n_elements(data_alt) emis_height_max = data_alt[n_alt-1] ; get wrfheight, wrflevs and dz (unit in meter) mean_height, emis_height_max, ph, phb, hgt, eta_dim, wrfheight, wrflevs, dz ;create empty WRF emissions netcdf files hrs = strcompress(string(indgen(24)),/remove_all) hrs[0:9] = '0'+hrs[0:9] ;timestring = yyyy+'-'+mm+'-'+dd+'_'+hrs+':00:00' ;file_out = 'wrfaircraftchemi_'+domain+'_'+timestring ;timestring = 'hour_'+hrs hrs = strcompress(string(indgen(24)),/remove_all) hrs[0:9] = '0'+hrs[0:9] timestring = yyyy+'-'+mm+'-'+dd+'_'+hrs+':00:00' file_out = 'wrfaircraftchemi_'+domain+'_'+timestring ncfiles = output_dir + file_out print, '' print, 'Creating the following files in ', output_dir for i = minhour, maxhour do begin $ create_netcdffile, ncfiles[i], timestring[i], wrflevs, wrf_ix, wrf_jx, dx*1e3, dx*1e3, lat1, lon1, truelat1, truelat2, moad_cen_lat, eta_dim, wrfheight, map_proj print,i+1,' ', file_out[i] endfor ; read aircraft emissions data ; spatial assignment aircraft -> WRFgrid (simple match by looking for the ; aircraft grid box nearest to a WRF pixel) openw,1,'data_WRFassign_dat' printf,1,wrf_ix, wrf_jx, wrflevs d=-99 dist = fltarr(wrf_ix, wrf_jx)-999. ; distance in meters between aircraft cell center and WRF cell center ; find the nearest horizontal data point index_lat = intarr(wrf_ix, wrf_jx) index_lon = intarr(wrf_ix, wrf_jx) for ix=0,wrf_ix-1 do begin for jx=0,wrf_jx-1 do begin klat = where (abs(wrf_lat[ix,jx]-data.lat) eq min(abs(wrf_lat[ix,jx]-data.lat))) klon = where (abs(wrf_lon[ix,jx]-data.lon) eq min(abs(wrf_lon[ix,jx]-data.lon))) distance,data.lat[klat[0]],data.lon[klon[0]],wrf_lat[ix,jx], wrf_lon[ix,jx], d dist[ix,jx] = d index_lat[ix,jx] = klat[0] index_lon[ix,jx] = klon[0] printf, 1,ix, jx, data.lat[klat[0]],data.lon[klon[0]],wrf_lat[ix,jx], wrf_lon[ix,jx], d*1e-3, format='(2I10, 5F10.3)' endfor endfor ; find the nearest vertical data point index_alt = intarr(wrflevs) ; data_alt is data.alt*1.e3 d_to_data_alt = data_alt for k=0,wrflevs-1 do begin for i=0,n_alt-1 do begin d_to_data_alt[i] = abs( wrfheight[k] - data_alt[i] ) endfor index_alt[k] = where ( d_to_data_alt eq min(d_to_data_alt) ) printf, 1, k, index_alt[k], format='(2I10)' endfor close,1 ;convert emission unit from molecules/cm3/sec to mole/km2/hour ;note: dz in meter arg_number = 6.022e23 ;;factor = 1.e5 * 1.e5 * 1.e2 *3600./arg_number factor = 1e5 * 1e5 *3600./arg_number for hrloop = minhour, maxhour do begin print, wrf_ix, wrf_jx, wrflevs dummy='' no_emis = fltarr(wrf_ix, wrf_jx, wrflevs) co_emis = fltarr(wrf_ix, wrf_jx, wrflevs) ch4_emis = fltarr(wrf_ix, wrf_jx, wrflevs) so2_emis = fltarr(wrf_ix, wrf_jx, wrflevs) for i = 0, wrf_ix-1 do begin for j = 0, wrf_jx-1 do begin for k = 0, wrflevs-1 do begin no_emis [i,j,k] = data.e_no [index_lon[i,j], index_lat[i,j], index_alt[k]] *1e5 * factor ;* dz[k] co_emis [i,j,k] = data.e_co [index_lon[i,j], index_lat[i,j], index_alt[k]] *1e5 * factor ;* dz[k] ch4_emis[i,j,k] = data.e_ch4[index_lon[i,j], index_lat[i,j], index_alt[k]] *1e5 * factor ;* dz[k] so2_emis[i,j,k] = data.e_so2[index_lon[i,j], index_lat[i,j], index_alt[k]] *1e5 * factor ;* dz[k] endfor ;plot, data.e_no [index_lon[i,j], index_lat[i,j],*] *factor * 1e-5, data.alt ;oplot,no_emis [i,j,*], wrfheight*1e-3, color=220 ;read,dummy endfor endfor ; set lowest km to zero (ac emissions included in anthropogenic ; emissions) ilev=where(wrfheight le maxheight) if ilev[0] ne -1 then begin no_emis[*,*,ilev] = 0. co_emis[*,*,ilev] = 0. ch4_emis[*,*,ilev] = 0. so2_emis[*,*,ilev] = 0. endif ; write to netcdf files ncid = ncdf_open(ncfiles[hrloop],/write) print,'Writing to ',ncfiles[hrloop] ncdf_varput,ncid,'EAC_NO', no_emis ncdf_varput,ncid,'EAC_CO', co_emis ncdf_varput,ncid,'EAC_CH4', ch4_emis ncdf_varput,ncid,'EAC_SO2', so2_emis ncdf_close,ncid endfor ; end of hour loop ;spawn, 'ncrcat wrfaircraftchemi_hour_0* wrfaircraftchemi_hour_10 wrfaircraftchemi_hour_11 wrfaircraftchemi_00z_d01' ;spawn, 'ncrcat wrfaircraftchemi_hour_1[23456789] wrfaircraftchemi_hour_2* wrfaircraftchemi_12z_d01' ;spawn, 'rm wrfaircraftchemi_hour_*' end ./PaxHeaders.31608/create_netcdffile.pro0000644000000000000000000000012712427253142015127 xustar0027 mtime=1415403106.664809 30 atime=1517528749.956060854 30 ctime=1415403106.664711201 create_netcdffile.pro0000755001151300017500000000541712427253142015522 0ustar00barthmncar00000000000000PRO create_netcdffile, emisfile, timestring, zdim, WE, SN, dx, dy, cenlat, cenlon, truelat1, truelat2, moad_cen_lat, eta_dim, z, mapproj ncid = ncdf_create(emisfile, /clobber) ; create dimensions timeid = ncdf_dimdef(ncid, 'Time',/UNLIMITED) weid = ncdf_dimdef(ncid, 'west_east',WE) snid = ncdf_dimdef(ncid, 'south_north',SN) btid = ncdf_dimdef(ncid, 'bottom_top',eta_dim) stagid = ncdf_dimdef(ncid,'ac_emissions_zdim_stag',zdim) datestrlenid = ncdf_dimdef(ncid, 'DateStrLen',19) ; write global attributes ncdf_attput,ncid,/GLOBAL,/long,'WEST-EAST_GRID_DIMENSION',WE+1 ncdf_attput,ncid,/GLOBAL,/long,'SOUTH-NORTH_GRID_DIMENSION',SN+1 ncdf_attput,ncid,/GLOBAL,/long,'BOTTOM-TOP_GRID_DIMENSION',eta_dim ncdf_attput,ncid,/GLOBAL,/char, 'TITLE','Created by Jeff Lee' ncdf_attput,ncid,/GLOBAL,/float,'DX',dx ncdf_attput,ncid,/GLOBAL,/float,'DY',dy ncdf_attput,ncid,/GLOBAL,/float,'HEIGHT',z ncdf_attput,ncid,/GLOBAL,/float,'CEN_LAT',cenlat ncdf_attput,ncid,/GLOBAL,/float,'CEN_LON',cenlon ncdf_attput,ncid,/GLOBAL,/float,'TRUELAT1',truelat1 ncdf_attput,ncid,/GLOBAL,/float,'TRUELAT2',truelat2 ncdf_attput,ncid,/GLOBAL,/float,'MOAD_CEN_LAT',moad_cen_lat ncdf_attput,ncid,/GLOBAL,/long,'MAP_PROJ',mapproj ; define varaiable fields tvarid = ncdf_vardef(ncid,'Times',[datestrlenid, timeid ],/char) varid = ncdf_vardef(ncid, 'EAC_NO', [weid,snid,stagid, timeid], /FLOAT) ncdf_attput,ncid,'EAC_NO', 'FieldType', /long, 104 ncdf_attput,ncid,'EAC_NO', 'MemoryOrder',/char,'XYZ' ncdf_attput,ncid,'EAC_NO', 'description',/char,'EMISSIONS' ncdf_attput,ncid,'EAC_NO', 'units', /char,'mole km-2 hr-1' ncdf_attput,ncid,'EAC_NO', 'stagger', /char,'Z' varid = ncdf_vardef(ncid, 'EAC_CO', [weid,snid,stagid, timeid], /FLOAT) ncdf_attput,ncid,'EAC_CO', 'FieldType', /long, 104 ncdf_attput,ncid,'EAC_CO', 'MemoryOrder',/char,'XYZ' ncdf_attput,ncid,'EAC_CO', 'description',/char,'EMISSIONS' ncdf_attput,ncid,'EAC_CO', 'units', /char,'mole km-2 hr-1' ncdf_attput,ncid,'EAC_CO', 'stagger', /char,'Z' varid = ncdf_vardef(ncid, 'EAC_CH4', [weid,snid,stagid, timeid], /FLOAT) ncdf_attput,ncid,'EAC_CH4', 'FieldType', /long, 104 ncdf_attput,ncid,'EAC_CH4', 'MemoryOrder',/char,'XYZ' ncdf_attput,ncid,'EAC_CH4', 'description',/char,'EMISSIONS' ncdf_attput,ncid,'EAC_CH4', 'units', /char,'mole km-2 hr-1' ncdf_attput,ncid,'EAC_CH4', 'stagger', /char,'Z' varid = ncdf_vardef(ncid, 'EAC_SO2', [weid,snid,stagid, timeid], /FLOAT) ncdf_attput,ncid,'EAC_SO2', 'FieldType', /long, 104 ncdf_attput,ncid,'EAC_SO2', 'MemoryOrder',/char,'XYZ' ncdf_attput,ncid,'EAC_SO2', 'description',/char,'EMISSIONS' ncdf_attput,ncid,'EAC_SO2', 'units', /char,'mole km-2 hr-1' ncdf_attput,ncid,'EAC_SO2', 'stagger', /char,'Z' ncdf_control,ncid, /endef ; write fields ncdf_varput,ncid,'Times',timestring ncdf_close,ncid end ./PaxHeaders.31608/ijll_lc.pro0000644000000000000000000000012612427253142013110 xustar0027 mtime=1415403106.669957 29 atime=1517528749.97430302 30 ctime=1415403106.669870411 ijll_lc.pro0000755001151300017500000000425012427253142013476 0ustar00barthmncar00000000000000PRO ijll_lc, i, j, lat, lon, truelat1, truelat2, hemi, stdlon,lat1, lon1, knowni, knownj, dx ; Subroutine to convert from the (i,j) cartesian coordinate to the ; geographical latitude and longitude for a Lambert Conformal projection. ; taken from WRF WPS geogrid rad_per_deg = !pi/180. deg_per_rad = 180/!pi rebydx = 6370./dx chi1 = (90. -hemi*truelat1)*rad_per_deg chi2 = (90. -hemi*truelat2)*rad_per_deg ; See if we are in the southern hemispere and flip the indices ; if we are. inew = hemi * i jnew = hemi * j ; Compute radius**2 to i/j location reflon = stdlon + 90. ala1 = lat1 * rad_per_deg alo1 = (lon1 - reflon) * rad_per_deg scale_top = 1. + hemi * SIN(truelat1 * rad_per_deg) deltalon1 = lon1 - stdlon IF (deltalon1 gt +180.) then deltalon1 = deltalon1 - 360. IF (deltalon1 lt -180.) then deltalon1 = deltalon1 + 360. IF (ABS(truelat1-truelat2) gt 0.1) THEN begin cone = ALOG10(COS(truelat1*rad_per_deg)) - $ ALOG10(COS(truelat2*rad_per_deg)) cone = cone /(ALOG10(TAN((45.0 - ABS(truelat1)/2.0) * rad_per_deg)) - $ ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg))) endif else if (ABS(truelat1-truelat2) le 0.1) then begin cone = SIN(ABS(truelat1)*rad_per_deg ) ENDIF tl1r = truelat1 * rad_per_deg ctl1r = COS(tl1r) rsw = rebydx * ctl1r/cone * (TAN((90.*hemi-lat1)*rad_per_deg/2.) / (TAN((90.*hemi-truelat1)*rad_per_deg/2.)))^cone arg = cone*(deltalon1*rad_per_deg) polei = hemi*knowni - hemi * rsw * SIN(arg) polej = hemi*knownj + rsw * COS(arg) xx = inew - polei yy = polej - jnew r2 = (xx*xx + yy*yy) r = SQRT(r2)/rebydx ;Convert to lat/lon IF (r2 eq 0.) THEN begin lat = hemi * 90. lon = stdlon endif else if r2 ne 0 then begin lon = stdlon + deg_per_rad * ATAN(hemi*xx,yy)/cone ; lon = AMOD(lon+360., 360.) lon = lon+360 MOD 360. IF (chi1 eq chi2) THEN begin chi = 2.0*ATAN( ( r/TAN(chi1) )^(1./cone) * TAN(chi1*0.5) ) endif ELSE if chi1 ne chi2 then begin chi = 2.0*ATAN( (r*cone/SIN(chi1))^(1./cone) * TAN(chi1*0.5)) ENDIF lat = (90.0-chi*deg_per_rad)*hemi ENDIF IF (lon gt +180.) then lon = lon - 360. IF (lon lt -180.) then lon = lon + 360. end ./PaxHeaders.31608/llij_lc.pro0000644000000000000000000000012712427253142013111 xustar0027 mtime=1415403106.690199 30 atime=1517528749.975829808 30 ctime=1415403106.690037429 llij_lc.pro0000755001151300017500000000271112427253142013476 0ustar00barthmncar00000000000000PRO llij_lc, xlat,xlon,xi,yj, truelat1, truelat2, hemi, stdlon,lat1, lon1, knowni, knownj, dx rad_per_deg = !pi/180. deg_per_rad = 180/!pi rebydx = 6370./dx deltalon = xlon - stdlon IF (deltalon gt +180.) then deltalon = deltalon - 360. IF (deltalon lt -180.) then deltalon = deltalon + 360. deltalon1 = lon1 - stdlon IF (deltalon1 gt +180.) then deltalon1 = deltalon1 - 360. IF (deltalon1 lt -180.) then deltalon1 = deltalon1 + 360. tl1r = truelat1 * rad_per_deg ctl1r = COS(tl1r) IF (ABS(truelat1-truelat2) gt 0.1) THEN begin cone = ALOG10(COS(truelat1*rad_per_deg)) - $ ALOG10(COS(truelat2*rad_per_deg)) cone=cone/(ALOG10(TAN((45.0-ABS(truelat1)/2.0)*rad_per_deg))- $ ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg))) endif else begin cone = SIN(ABS(truelat1)*rad_per_deg ) ENDelse rm = rebydx * ctl1r/cone * $ (TAN((90.*hemi-xlat)*rad_per_deg/2.) / $ TAN((90.*hemi-truelat1)*rad_per_deg/2.))^cone arg = cone*(deltalon*rad_per_deg) arg1 = cone*(deltalon1*rad_per_deg) rsw = rebydx * ctl1r/cone * $ (TAN((90.*hemi-lat1)*rad_per_deg/2.) / $ TAN((90.*hemi-truelat1)*rad_per_deg/2.))^cone polei = hemi*knowni - hemi * rsw * SIN(arg1) polej = hemi*knownj + rsw * COS(arg1) xi = polei + hemi * rm * SIN(arg) yj = polej - rm * COS(arg) xi = hemi * xi yj = hemi * yj END ./PaxHeaders.31608/mapcf.pro0000644000000000000000000000012712427253142012567 xustar0027 mtime=1415403106.708593 30 atime=1517528749.993136227 30 ctime=1415403106.708520751 mapcf.pro0000755001151300017500000000174612427253142013163 0ustar00barthmncar00000000000000; adapted from fortran suroutine MAPCF by G. Frost ; computes the lat and lon from model indexes of ; a point; Lambert Conformal Projection ; Gabriele Pfister 2004 ; this one is needed for NEI emissions PRO MAPCF, xi, yj, xlat, xlon dxkm = 4. xlatc = 40. xlonc = -97. il = 1369. jl = 1045 clat1 = 45. clat2 = 33. A = 6370.997 pole = 90. count=0L frty5d = atan(1.) conv = 45./frty5d rlat1 = clat1/conv rlat2 = clat2/conv xn = alog(cos(rlat1)/cos(rlat2))/alog(tan(frty5d+.5*rlat2)/tan(frty5d+.5*rlat1)) psx = (pole-xlatc)/conv cell = a*sin(rlat1)/xn cell2 = (tan(psx/2.))/(tan(rlat1/2.)) R = cell*(cell2^xn) cntrj = (jl+1)/2. cntri = (il+1)/2. rxn = 1.0/xn x = (xi-cntri)*dxkm y = (yj-cntrj)*dxkm-R FLP = ATAN(-1*x/y) FLPP = (FLP/XN)*CONV+xlonc IF (FLPP LT -180) THEN FLPP = FLPP + 360. IF (FLPP GT 180) THEN FLPP = FLPP-360. xlon=FLPP zzzz = x*y+y*y rs = SQRT(X*X+y*Y) cell = (rs*xn)/(a*sin(rlat1)) cel1 = tan(rlat1/2.)*(cell^rxn) cel2 = atan(cel1) psx = 2*cel2*conv xlat = pole-psx END ./PaxHeaders.32448/README0000644000000000000000000000012413235104537011640 xustar0027 mtime=1517586783.470656 27 atime=1517586783.470484 30 ctime=1517586783.470598917 README0000644001151300017500000000254713235104537012231 0ustar00barthmncar00000000000000Aircraft Emissions Preprocessor Tool This IDL code reads the netcdf file baughcum_1999.aircraft.1x1.nc and interpolates emissions onto the WRF grid of interest. The "user modify" section is contained in lines 135-165. In this section, one defines where the wrfinput file and data file are located, the name of the domain (d01, d02, d03, etc), date (year, month, day), the time of day for the files (hourly output for 24 hours is the current set up), and the maximum height of the anthropogenic emissions. This maximum height for anthropogenic emissions is needed because the US EPA NEI contains aircraft emissions and therefore we set aircraft emissions in the output here to zero for heights below the maximum height. Files created are named 'wrfaircraftchemi_'+domain+'_'+yyyy+'-'+mm+'-'+dd+'_'+hrs+':00:00'. In other words, the same naming protocol is used as for wrfchemi files. The code is currently set up to create 24 files, one for each hour of day. If 2 12-hour files are desired, there are lines at the bottom of the script that are commented out to show you how to combine the individual files into 2 12-hour files using nco commands. That is, /usr/bin/ncrcat wrfaircraftchemi_hour_0* wrfaircraftchemi_hour_10 wrfaircraftchemi_hour_11 wrfaircraftchemi_00z_d02 /usr/bin/ncrcat wrfaircraftchemi_hour_1[23456789] wrfaircraftchemi_hour_2* wrfaircraftchemi_12z_d01