;************************************ ; derive.ncl ; derives various hao quantities. ;************************************* undef("get_mol_wt") function get_mol_wt(name:string) local wt begin wt = new(1,"float") if(name.eq."ARGON"); AR wt = 40. return(wt) end if if(name.eq."CH4") wt = 16. return(wt) end if if(name.eq."CARBON MONOXIDE") ; CO wt = 28. return(wt) end if if(name.eq."CARBON DIOXIDE") ; CO2 wt = 44. return(wt) end if if(name.eq."HYDROGEN");H wt = 1. return(wt) end if if(name.eq."MOLECULAR HYDROGEN"); H2 wt = 2. return(wt) end if if(name.eq."WATER VAPOR") ; H20 wt = 18. return(wt) end if if(name.eq."H2P") wt = 2.2 return(wt) end if if(name.eq."H2V1".or.name.eq."H2V2".or.name.eq."H2V3".or.name.eq."H2V4")then wt = 2.2 return(wt) end if if(name.eq."H3P") wt = 3. return(wt) end if if(name.eq."HELIUM"); HE wt = 4. return(wt) end if if(name.eq."HO2") wt = 33. return(wt) end if if(name.eq."HYDROGEN GROUP (OH+HO2+H)") ; HOX wt = 17. return(wt) end if if(name.eq."HP") wt = 1. return(wt) end if if(name.eq."N2") wt = 28. return(wt) end if if(name.eq."N2D");N2D wt = 14. return(wt) end if if(name.eq."N4S") ; NS4 wt = 14. return(wt) end if if(name.eq."TOTAL SODIUM"); NAT wt = 23. return(wt) end if if(name.eq."NITRIC OXIDE"); NO wt = 30. return(wt) end if if(name.eq."NITROGEN DIOXIDE"); NO2 wt = 46. return(wt) end if if(name.eq."NO GROUP (NO+NO2)") ;NOZ wt = 30. return(wt) end if if(name.eq."ATOMIC OXYGEN");O1 wt = 16. return(wt) end if if(name.eq."MOLECULAR OXYGEN"); O2 wt = 32. return(wt) end if if(name.eq."O21D"); O21D wt = 32. return(wt) end if if(name.eq."OZONE"); ozone wt = 48. return(wt) end if if(name.eq."HYDROXYL"); OH wt = 17. return(wt) end if if(name.eq."OX GROUP (O+O3)") ; OX wt = 16. return(wt) end if print("Variable "+name+" does not have a molecular weight listed.") print("It is not a density variable. Can not convert units") ; wt = -999. return(wt) end ;*************************************************** undef("den_convert_units") function den_convert_units(x:numeric,in:file,units:string) local barm,pkt,ntime,nlev,nlat,nlon,dims,comp,w,zp,zpbot,zptop,O1,O2,boltz,dz,TN begin if(.not.isDensity(x@name))then ; this function searches by short_name ; print("not a density variable, skipping conversion") return(x) end if if(dimsizes(dimsizes(x)).ne.4)then return(x) end if boltz = 1.3805e-16 p0 = 5.e-4 dims = dimsizes(x) nlev = dims(1) ; calculate barm O2 = in->O2 O1 = in->O1 ; barm = new( (/ntime,nlev,nlat,nlon/),typeof(x) ) barm = x N2 = (1-O2-O1) ; barm is the mean molecular weight barm = 1./(O2/32.+O1/16.+N2/28.) ; calculate pkt ; pkt = new( (/ntime,nlev,nlat,nlon/),typeof(x) ) pkt = x zp = x&lev if(isfilevar(in,"TN"))then TN = in->TN else print("TN is required for density conversions") print("TN not on file") return(x) end if do k = 0,nlev-1 tmp = exp(-zp(k)) pkt(:,k,:,:) = p0*tmp/(boltz*TN(:,k,:,:)) end do ; get molecular weight wt = get_mol_wt(x@long_name) if(ismissing(wt))then return(x) ; not a density variable. don't convert end if ; begin conversions if(units.eq."CM3")then print("converting "+x@long_name+" to CM3 units") x = (/ x*pkt*barm/wt /) x@units = "cm3" end if if(units.eq."CM3-MR")then print("converting "+x@long_name+" to CM3-MR units") x = (/ x*barm/wt /) x@units = "cm3-MR" end if if(units.eq."GM/CM3")then print("converting "+x@long_name+" to GM/CM3 units") x = (/ x*pkt*barm*1.66e-24/) x@units = "GM/cm3" end if if(units.eq."MMR")then print("no conversion required") end if return(x) end ;************************************************ undef("derive") function derive(name:string,in:file,units:string) local rho,O2,O1,N2,x,CO,CO2,CH4,ctot ; available derivations: CTOT,N2,N2/O1,O1/N2,O1/O2+N2,O1/O2,RHO begin if(name.eq."CTOT")then ; CO+CO2+CH4 if(isfilevar(in,"CO"))then if(isfilevar(in,"CO2"))then if(isfilevar(in,"CH4"))then ; read in variables and convert to full zp CO = in->CO CO = halfzp("CO",CO) CO@short_name = "CO" CO = den_convert_units(CO,in,units) CO2 = in->CO2 CO2 = halfzp("CO2",CO2) CO2@short_name = "CO2" CO2 = den_convert_units(CO2,in,units) CH4 = in->CH4 CH4 = halfzp("CH4",CH4) CH4@short_name = "CH4" CH4 = den_convert_units(CH4,in,units) ; derive main variable ctot = CO ctot = CO+CO2+CH4 ctot@long_name = "Total Carbon" ; fix bad units ctot = coordcheck(ctot) else print("CH4 not on file, can not derive CTOT") return end if else print("CO2 not on file, can not derive CTOT") return end if else print("CO not on file, can not derive CTOT") return end if return(ctot) end if ;******************************** if(name.eq."RHO")then ; O2+O1+N2 where N2=1-O2-O1 if(isfilevar(in,"O2"))then if(isfilevar(in,"O1"))then ; read in variables and convert to full zp O2 = in->O2 O2 = halfzp("O2",O2) O1 = in->O1 O1 = halfzp("O1",O1) ; derive sub-variables(this must be done with 02 & 01 in MMR N2 = O2 ; trick to retain coordinate info N2 = 1-O2-O1 N2@short_name = "N2" N2 = den_convert_units(N2,in,units) ; convert before calculation O2@short_name = "O2" O2 = den_convert_units(O2,in,units) O1@short_name = "O1" O1 = den_convert_units(O1,in,units) ; derive main variable rho = O2 ; trick to retain coordinate info rho = O2+O1+N2 rho@long_name = "RHO" rho@short_name = "RHO" ; fix bad units rho = coordcheck(rho) else print("O1 not on file, can not derive RHO") return end if else print("O2 not on file, can not derive RHO") return end if return(rho) end if ;*********************** if(name.eq."N2")then ; N2 = 1-O2-O1 if(isfilevar(in,"O2"))then if(isfilevar(in,"O1"))then ; read in variables and convert to full zp O2 = in->O2 O2 = halfzp("O2",O2) O1 = in->O1 O1 = halfzp("O1",O1) ; derive sub-variables(this must be done with 02 & 01 in MMR N2 = O2 ; trick to retain coordinate info N2 = 1-O2-O1 N2@short_name = "N2" N2 = den_convert_units(N2,in,units) ; fix bad units N2 = coordcheck(N2) else print("O1 not on file, can not derive N2") return end if else print("O2 not on file, can not derive N2") return end if return(N2) end if ;************************************** if(name.eq."O1/N2")then ; O1/N2 where N2=1-O2-O1 ; O1 is converted to units prior to ratio if(isfilevar(in,"O2"))then if(isfilevar(in,"O1"))then ; read in variables and convert to full zp O2 = in->O2 O2 = halfzp("O2",O2) O1 = in->O1 O1 = halfzp("O1",O1) ; derive sub-variables N2 = O2 N2 = 1-O2-O1 N2@short_name = "N2" N2 = den_convert_units(N2,in,units) O1@short_name = "O1" O1 = den_convert_units(O1,in,units) ; derive main variable x = O1 x = O1/N2 x@long_name = "ratio O1/N2" ; fix bad units x = coordcheck(x) else print("O1 not on file, can not derive O1/N2") return end if else print("O2 not on file, can not derive O1/N2") return end if return(x) delete(x) end if ;************************************* if(name.eq."N2/O1")then ; O1/N2 where N2=1-O2-O1 if(isfilevar(in,"O2"))then if(isfilevar(in,"O1"))then ; read in variables and convert to full zp O2 = in->O2 O2 = halfzp("O2",O2) O1 = in->O1 O1 = halfzp("O1",O1) ; derive sub-variables N2 = O2 N2 = 1-O2-O1 N2@short_name = "N2" N2 = den_convert_units(N2,in,units) O2@short_name = "O2" O2 = den_convert_units(O2,in,units) ; derive main variable x = O1 x = N2/O1 x@long_name = "ratio N2/O1" ; fix bad units x = coordcheck(x) else print("O1 not on file, can not derive N2/O1") return end if else print("O2 not on file, can not derive N2/O1") return end if return(x) delete(x) end if ;********************************* if(name.eq."O1/O2+N2")then ; O1/N2 where N2=1-O2-O1 if(isfilevar(in,"O2"))then if(isfilevar(in,"O1"))then ; read in variables and convert to full zp O2 = in->O2 O2 = halfzp("O2",O2) O1 = in->O1 O1 = halfzp("O1",O1) ; derive sub-variables N2 = O2 N2 = 1-O2-O1 N2@short_name = "N2" N2 = den_convert_units(N2,in,units) O2@short_name = "O2" O2 = den_convert_units(O2,in,units) O1@short_name = "O1" O1 = den_convert_units(O1,in,units) ; derive main variable x = O1 x = O1/(O2+N2) x@long_name = "ratio O1/O2+N2" ; fix bad units x = coordcheck(x) else print("O1 not on file, can not derive O1/O2+N2") return end if else print("O2 not on file, can not derive O1/O2+N2") return end if return(x) delete(x) end if ;***************************** if(name.eq."O1/O2")then if(isfilevar(in,"O2"))then if(isfilevar(in,"O1"))then ; read in variables and convert to full zp O2 = in->O2 O2 = halfzp("O2",O2) O2@short_name = "O2" O2 = den_convert_units(O2,in,units) O1 = in->O1 O1 = halfzp("O1",O1) O1@short_name = "O1" O1 = den_convert_units(O1,in,units) ; derive main variable x = O1 x = O1/O2 x@long_name = "ratio O1/O2" ; fix bad units x = coordcheck(x) else print("O1 not on file, can not derive O1/O2") return end if else print("O2 not on file, can not derive O1/O2") return end if return(x) delete(x) end if end ;****************************************