####WRF CHEM output into PM25 and PM10####

1. To get various variables from WRF OUTPUT following python codes and references were used #####to get WRF output NETCDF into python##### The below code import netcdf file into python and subset the data based on Coimbatore domain of

from netCDF4 import Dataset
import numpy as np
wrfoutput = 'wrfout_d01_2014-06-05_05:00:00_D03'
fh = Dataset(wrfoutput, mode='r')
#to view all the variables in imported wrf output netcdf
vars = fh.variables
for v in vars:
	print v

1. To subset the output into Coimbatore domain of Left-down corner (76.7900,10.8614), left up corner (76.7900, 11.2300), right down corner(76.1600,10.8614) right up corner (76.1600,11.2300). ``` lons = fh.variables[‘XLONG’][:] lats = fh.variables[‘XLAT’][:] #to print variable names with its shape vars = fh.variables for v in vars: lons = fh.variables[v][:] d=lons.shape print str(v)+’:‘+str(d)

#to reshape the numpy array from one dimension to another lats.shape (1,72,75) x = lats.reshape((72,75)) #to get the index of array values y=np.where(x==76.1600)

1. The complete code of script to convert the wrf output into pandas dataframe with subset into Coimbatore domain

from netCDF4 import Dataset import numpy as np wrfoutput = ‘wrfoutd012014-06-0505:00:00D03’ fh = Dataset(wrfoutput, mode=‘r’) lons = fh.variables[‘XLONG’][:] lats = fh.variables[‘XLAT’][:] lons.shape lonsX0=lons.reshape((105,96)) lonsX=lonsX0.reshape((10080)) cbeX=lonsX[333:1666] #the above value is from the all lat and long #values was converted into pandas dataframe and then csv # #and viewed for location of bounding box lat long #(76.16,10.86) and (76.79, 11.23), it is in located in the #index beween 333:1666 latsX0=lats.reshape((105,96)) latsX=latsX0.reshape((10080)) cbeY=latsX[333:1666]

####Temperature at 2 m height###

1. to get the temperrature at 2m variable from wrf out put based on [this](http://www.meteo.unican.es/wiki/cordexwrf/OutputVariables) regarding the variable name in netcdf. To get the variable in 

Temp =  fh.variables['T2'][:]
Temp0=Temp.reshape((105,96))
Temp1=Temp0.reshape((10080))
cbeT = Temp1[333:1666]

####Relative Humidity at 2 m Height####

Based on [this](http://mailman.ucar.edu/pipermail/wrf-users/2012/002546.html)

psfc = fh.variables[‘PSFC’][:] psfc0=psfc.reshape((105,96)) psfc1=psfc.reshape((10080)) cbePSFC = psfc1[333:1666] q2=fh.variables[‘Q2’][:] q20=q2.reshape((105,96)) q21=q20.reshape((10080)) cbeQ2 = q21[333:1666] pq0 = 379.90516 a2 = 17.2693882 a3 = 273.16 a4 = 35.86 #The formula [f_rh2 = q2 / ( (pq0 / psfc) * exp(a2 * (t2 - #a3) / (t2 - a4)) )] d=cbeQ2/((379.90516/cbePSFC)(np.exp(17.2693882(cbeT-273.16/cbeT-35.86)))

This was diffcult to do and found a easier method using ncl script to get not only RH but also Temp, WS, WD.

####NCL script to convert the wrf output to get the variables, TEMP, RH, WS, WD####
Based on [this](http://foehn.colorado.edu/wrfout_to_cf/variables.html)
1. to run the script used following command

```ncl 'file_in="wrfout_d01_2014-06-05_00:00:00_D03"' 'file_out="wrfpost.nc"' wrfout_to_cf.ncl```

from netCDF4 import Dataset
import numpy as np
wrfoutput = 'wrfpost.nc'
fh = Dataset(wrfoutput, mode='r')
temp = fh.variables['T_2m'][:]
Temp0=temp.reshape((105,96))
Temp1=Temp0.reshape((10080))
cbeT = Temp1[333:1666]

hum = fh.variables['rh_2m'][:]
hum0=hum.reshape((105,96))
hum1=hum0.reshape((10080))
cbeRH = hum1[333:1666]

ws = fh.variables['ws_10m'][:]
ws0=ws.reshape((105,96))
ws1=ws0.reshape((10080))
cbeWS = ws1[333:1666]

wd = fh.variables['wd_10m'][:]
wd0=wd.reshape((105,96))
wd1=wd0.reshape((10080))
cbeWD = wd1[333:1666]

####PM25 and PM10 from wrf output####
1. Based on this [paper](http://www.geosci-model-dev.net/7/1621/2014/gmd-7-1621-2014.pdf)
for PM2.5

PM25= pd[p25+1.375S+BC1+BC2+1.8(OC1+OC2)+D1+0.286D2+SS1+0.942SS2]

PM10=pd[p25+1.375S+BC1+BC2+1.8(OC1+OC2)+D1+D2+D3+0.87D4+SS1+SS2+SS3]

wrfoutput = ‘wrfoutd012014-06-0505:00:00D03’ fh = Dataset(wrfoutput, mode=‘r’) for v in vars: lons = fh.variables[v][:] d=lons.shape print str(v)+’:‘+str(d)

BC1 = fh.variables[‘BC1’][:] BC1=BC1[:,0] BC1r=BC1.reshape((105,96)) BC1d=BC1r.reshape((10080)) cbeBC1 = BC1d[333:1666]

BC2 = fh.variables[‘BC2’][:] BC2=BC2[:,0] BC2r=BC2.reshape((105,96)) BC2d=BC2r.reshape((10080)) cbeBC2 = BC2d[333:1666]

OC1 = fh.variables[‘OC1’][:] OC1=OC1[:,0] OC1r=OC1.reshape((105,96)) OC1d=OC1r.reshape((10080)) cbeOC1 = OC1d[333:1666]

OC2 = fh.variables[‘OC2’][:] OC2=OC2[:,0] OC2r=OC2.reshape((105,96)) OC2d=OC2r.reshape((10080)) cbeOC2 = OC2d[333:1666]

DUST1 = fh.variables[‘DUST1’][:] DUST1=DUST_1[:,0] DUST1r=DUST1.reshape((105,96)) DUST1d=DUST1r.reshape((10080)) cbeD1 = DUST1d[333:1666]

DUST2 = fh.variables[‘DUST2’][:] DUST2=DUST_2[:,0] DUST2r=DUST2.reshape((105,96)) DUST2d=DUST2r.reshape((10080)) cbeD2 = DUST2d[333:1666]

DUST3 = fh.variables[‘DUST3’][:] DUST3=DUST_3[:,0] DUST3r=DUST3.reshape((105,96)) DUST3d=DUST3r.reshape((10080)) cbeD3 = DUST3d[333:1666]

DUST4 = fh.variables[‘DUST4’][:] DUST4=DUST_4[:,0] DUST4r=DUST4.reshape((105,96)) DUST4d=DUST4r.reshape((10080)) cbeD4 = DUST4d[333:1666]

SEAS1 = fh.variables[‘SEAS1’][:] SEAS1=SEAS_1[:,0] SEAS1r=SEAS1.reshape((105,96)) SEAS1d=SEAS1r.reshape((10080)) cbeSS1 = SEAS1d[333:1666]

SEAS2 = fh.variables[‘SEAS2’][:] SEAS2=SEAS_2[:,0] SEAS2r=SEAS2.reshape((105,96)) SEAS2d=SEAS2r.reshape((10080)) cbeSS2 = SEAS2d[333:1666]

SEAS3 = fh.variables[‘SEAS3’][:] SEAS3=SEAS_3[:,0] SEAS3r=SEAS3.reshape((105,96)) SEAS3d=SEAS3r.reshape((10080)) cbeSS3 = SEAS3d[333:1666]

Temp = fh.variables[‘T2’][:] #Temp=Temp[:,0] Temp0=Temp.reshape((105,96)) Temp1=Temp0.reshape((10080)) cbeT = Temp1[333:1666]

PSFC = fh.variables[‘PSFC’][:] PSFC=PSFC.reshape((105,96)) PSFC=PSFC.reshape((10080)) cbePSFC = PSFC[333:1666]

a=287.05*cbeT pd= cbePSFC/a

P25 = fh.variables[‘P25’][:] P25=P25[:,0] P25=P25.reshape((105,96)) P25=P25.reshape((10080)) cbeP25 = P25[333:1666]

sulf = fh.variables[‘sulf’][:] sulf=sulf[:,0] sulf=sulf.reshape((105,96)) sulf=sulf.reshape((10080)) cbeS = sulf[333:1666]

PM25= pd(cbeP25+1.375cbeS+cbeBC1+cbeBC2+1.8(cbeOC1+cbeOC2)+beD1+0.286cbeD2+cbeSS1+0.0942*cbeSS2)

PM10= pd(cbeP25+(1.375cbeS)+cbeBC1+cbeBC2+(1.8(cbeOC1+cbeOC2))+cbeD1+cbeD2+cbeD3+(0.87cbeD4)+cbeSS1+cbeSS2+cbeSS3)

P10 = fh.variables[‘P10’][:] P10=P10[:,0] P10=P10.reshape((105,96)) P10=P10.reshape((10080)) cbeP10 = P10[333:1666]

data1=pd.DataFrame(cbeX) data1[‘cbeY’]=pd.DataFrame(cbeY) data1[‘cbeT’]=pd.DataFrame(cbeT) data1[‘cbeRH’]=pd.DataFrame(cbeRH) data1[‘cbeWS’]=pd.DataFrame(cbeWS) data1[‘cbeWD’]=pd.DataFrame(cbeWD) data1[‘cbeP25’]=pd.DataFrame(PM25) data1[‘cbeP10’]=pd.DataFrame(PM10)

WIND DIRECTION based on this

https://www.eol.ucar.edu/content/wind-direction-quick-reference

Dirmet = atan2(-Umet,-Vmet) * DperR = 270 - ( atan2(Vmet,Umet) * DperR ) Spd = sqrt(Umet2 + Vmet2) Umet = -Spd * sin(Dirmet * RperD) Vmet = -Spd * cos(Dirmet * RperD)

wind earth rotation

Uearth = Ucosalpha + Vsinalpha, Vearth = Vcosalpha - Usinalpha.