PRO tgcm_sf_nocool ; read 5 minutely sampled TGCM netCDF output ;dtr = !pi/180. ;Missing value = 1.e36 missing = 1.e36 ; p0 is not a parameter, as it must be changed for mtgcm (see getflds) ; (for tgcm, p0=5.e-4, for mtgcm p0=1.2e-3) ; p0 = 5.e-4 ;microbar grav = 870. ; cm/s2 ;k0 is Boltzmann constant k0 = 1.38e-16 ; gas constant gask = 8.314e7 ;R ;Ri is Earth's radius at 0 km altitude in cm Ri = 6370*1.e5 dlat = 5.*!pi/180 dlon = 5.*!pi/180 ;file 11111111111111111111111111------------------------------------------------------------------- ;files = file_search('$LINUX/SolarFlare/sfism_all_yr200331*.nc') files = file_search('/hao/aim3/huangys/tiegcm1.94_data/SolarFlare/nc2003301/sfism_preflare_2003301.nc') nfiles = n_elements(files) print, 'number of files = ',nfiles mtime = intarr(3,288*nfiles) mut = fltarr(288*nfiles) mhpower = fltarr(288*nfiles) msigma_ped = fltarr(72,36,29,288*nfiles);S/m meden = fltarr(72, 36, 29, 288*nfiles) ;cm-3 mden = fltarr(72,36,29,288*nfiles);ilev mqjoule = fltarr(72,36,29,288*nfiles) ; (lon,lat,lev,time) mqjoule_integ = fltarr(72,36,288*nfiles); (lon,lat,time) msheating = fltarr(72,36,29,288*nfiles) ;mo1 = fltarr(72,36,29,288*nfiles) ;mo2 = fltarr(72,36,29,288*nfiles) mtn = fltarr(72,36,29,288*nfiles) ;mscht = fltarr(72,36,29,288*nfiles) mzg = fltarr(72,36,29,288*nfiles);ilev mz = fltarr(72,36,29,288*nfiles);ilev mtec = fltarr(72,36,288*nfiles) mqop2pa = fltarr(72,36,29,288*nfiles) ; production rate (cm-3s-1) mqop2da = fltarr(72,36,29,288*nfiles) ; (lon,lat,lev,time) mqo2p = fltarr(72,36,29,288*nfiles) ; (lon,lat,lev,time) mqn2p = fltarr(72,36,29,288*nfiles) ; (lon,lat,lev,time) mqnp = fltarr(72,36,29,288*nfiles) ; (lon,lat,lev,time) mqopa = fltarr(72,36,29,288*nfiles) ; (lon,lat,lev,time) mqnop = fltarr(72,36,29,288*nfiles) ; (lon,lat,lev,time) mco2cool = fltarr(72,36,29,288*nfiles) mnocool = fltarr(72,36,29,288*nfiles) for i = 0L,nfiles-1 do begin print, 'i=',i filein = files(i) ;; open and read netCDF file print, 'Reading ',filein id=ncdf_open(filein) ncdf_varget,id,'mtime',var mtime(*,i*288:(i*288+287)) = var ncdf_varget,id,'lat',lat ;(36),[-87.5,87.5] ncdf_varget,id,'lon',lon ;(72),[-180., 175.] ncdf_varget,id,'lev',lev ; ncdf_varget,id,'ilev',ilev ;Z,ZG,DEN ncdf_varget,id,ncdf_varid(id,'CO2_COOL'),vco2 ;erg/g/s mco2cool(*,*,*,i*288:(i*288+287)) = vco2 ncdf_varget,id,ncdf_varid(id,'NO_COOL'),vno ;erg/g/s mnocool(*,*,*,i*288:(i*288+287)) = vno ncdf_varget,id,ncdf_varid(id,'ut'),vut ;universal time (hours) mut(i*288:(i*288+287)) = vut ncdf_varget,id,ncdf_varid(id,'hpower'),vhpower ;hemispheric power, GW mhpower(i*288:(i*288+287)) = vhpower ncdf_varget,id,ncdf_varid(id,'SIGMA_PED'),vped ;pederson conductivity msigma_ped(*,*,*,i*288:(i*288+287)) = vped ncdf_varget,id,ncdf_varid(id,'NE'),veden ;electron density, /cm3 meden(*,*,*,i*288:(i*288+287)) = veden ncdf_varget,id,ncdf_varid(id,'DEN'),vden ;total density, g/cm3 mden(*,*,*,i*288:(i*288+287)) = vden ncdf_varget,id,ncdf_varid(id,'QJOULE'),vqj mqjoule(*,*,*,i*288:(i*288+287)) = vqj ncdf_varget,id,ncdf_varid(id,'QJOULE_INTEG'),vqj_int ;erg/cm2/s = mW/m2 mqjoule_integ(*,*,i*288:(i*288+287)) = vqj_int ncdf_varget,id,ncdf_varid(id,'SHEATING'),vheat ;solar heating,erg/g/s msheating(*,*,*,i*288:(i*288+287)) = vheat ;ncdf_varget,id,ncdf_varid(id,'QTOTAL'),qtotal ; ncdf_varget,id,ncdf_varid(id,'O1'),vo1 ;mo1(*,*,*,i*288:(i*288+287)) = vo1 ; ncdf_varget,id,ncdf_varid(id,'O2'),vo2 ;mo2(*,*,*,i*288:(i*288+287)) = vo2 ncdf_varget,id,ncdf_varid(id,'TN'),vtn mtn(*,*,*,i*288:(i*288+287)) = vtn ncdf_varget,id,ncdf_varid(id,'grav'),grav ;cm/s ; ncdf_varget,id,ncdf_varid(id,'SCHT'),vsh ;scale height, km ;mscht(*,*,*,i*288:(i*288+287)) = vsh ncdf_varget,id,ncdf_varid(id,'ZG'),vzg mzg(*,*,*,i*288:(i*288+287)) = vzg ncdf_varget,id,ncdf_varid(id,'Z'),vz mz(*,*,*,i*288:(i*288+287)) = vz ncdf_varget,id,ncdf_varid(id,'TEC'),vtec mtec(*,*,i*288:(i*288+287)) = vtec ;1/cm2 ncdf_varget,id,ncdf_varid(id,'QOP2Pa'),vqop2pa mqop2pa(*,*,*,i*288:(i*288+287)) = vqop2pa ncdf_varget,id,ncdf_varid(id,'QOP2Da'),vqop2da mqop2da(*,*,*,i*288:(i*288+287)) = vqop2da ncdf_varget,id,ncdf_varid(id,'QO2P'),vqo2p mqo2p(*,*,*,i*288:(i*288+287)) = vqo2p ncdf_varget,id,ncdf_varid(id,'QN2P'),vqn2p mqn2p(*,*,*,i*288:(i*288+287)) = vqn2p ncdf_varget,id,ncdf_varid(id,'QNP'),vqnp mqnp(*,*,*,i*288:(i*288+287)) = vqnp ncdf_varget,id,ncdf_varid(id,'QOPa'),vqopa mqopa(*,*,*,i*288:(i*288+287)) = vqopa ncdf_varget,id,ncdf_varid(id,'QNOP'),vqnop mqnop(*,*,*,i*288:(i*288+287)) = vqnop print,'Reading Finished~~~~~~~~~~~~~~~' endfor mtime = transpose(mtime) ntimes = n_elements(mtime(*,0)) ;# of hours smtime = fltarr(ntimes) ;for i = 0, ntimes-1 do begin smtime(*) = mtime(*,0)+mtime(*,1)/24.+mtime(*,2)/1440. ;smtime(i) = mtime(i*8,0) ;daily average ;endfor nlat = n_elements(lat) nlon = n_elements(lon) nlev = n_elements(lev) fp = where(mco2cool ge missing) mco2cool(fp) = !Values.F_NAN ;NaN fp = where(mnocool ge missing) mnocool(fp) = !Values.F_NAN ;NaN fp = where(mhpower ge missing) mhpower(fp) = !Values.F_NAN ;NaN fp = where(msigma_ped ge missing) msigma_ped(fp) = !Values.F_NAN ;NaN fp = where(meden ge missing) meden(fp) = !Values.F_NAN ;NaN fp = where(mden ge missing) mden(fp) = !Values.F_NAN ;NaN fp = where(mqjoule ge missing) mqjoule(fp) = !Values.F_NAN ;NaN fp = where(mqjoule_integ ge missing) mqjoule_integ(fp) = !Values.F_NAN ;NaN fp = where(msheating ge missing) msheating(fp) = !Values.F_NAN ;NaN ;fp = where(mo1 ge missing) ;mo1(fp) = !Values.F_NAN ;NaN ;fp = where(mo2 ge missing) ;mo2(fp) = !Values.F_NAN ;NaN fp = where(mtn ge missing) mtn(fp) = !Values.F_NAN ;NaN ;fp = where(mscht ge missing) ;mscht(fp) = !Values.F_NAN ;NaN fp = where(mzg ge missing) mzg(fp) = !Values.F_NAN ;NaN fp = where(mz ge missing) mz(fp) = !Values.F_NAN ;NaN fp = where(mtec ge missing) mtec(fp) = !Values.F_NAN ;NaN fp = where(mqop2pa ge missing) mqop2pa(fp) = !Values.F_NAN ;NaN fp = where(mqop2da ge missing) mqop2da(fp) = !Values.F_NAN ;NaN fp = where(mqo2p ge missing) mqo2p(fp) = !Values.F_NAN ;NaN fp = where(mqn2p ge missing) mqn2p(fp) = !Values.F_NAN ;NaN fp = where(mqnp ge missing) mqnp(fp) = !Values.F_NAN ;NaN fp = where(mqopa ge missing) mqopa(fp) = !Values.F_NAN fp = where(mqnop ge missing) mqnop(fp) = !Values.F_NAN ;NaN mtpr = mqop2pa+mqop2da+mqo2p+mqn2p+mqnp+mqopa+mqnop ;total production rate ;----------------------------Finish of data loading ;nocooling @ 120km at different time ;GeometricHeight = 4.e7 ;cm GeometricHeight = 1.2e7 ;1.e7+findgen(30)*(550.-100.)/30*1.e5 ;cm nzg = n_elements(GeometricHeight) ; == 1 sample_time = 5./60./24.+301. tn = fltarr(nlon,nlat,nzg,ntimes) gtn = fltarr(ntimes,nzg) zg = fltarr(nlon,nlat,nzg,ntimes) gzg = fltarr(ntimes,nzg) den = fltarr(nlon,nlat,nzg,ntimes) gden = fltarr(ntimes,nzg) eden = fltarr(nlon,nlat,nzg,ntimes) geden = fltarr(ntimes,nzg) co2cool = fltarr(nlon,nlat,nzg,ntimes) gco2cool = fltarr(ntimes,nzg) nocool = fltarr(nlon,nlat,nzg,ntimes) gnocool = fltarr(ntimes,nzg) sheating=fltarr(nlon,nlat,nzg,ntimes) ; solar heating per volume sheating_pm=fltarr(nlon,nlat,nzg,ntimes) ; solar heating per mass gsheating=fltarr(ntimes,nzg) gmsheating=fltarr(ntimes,nzg) jheating=fltarr(nlon,nlat,nzg,ntimes) gjheating=fltarr(ntimes,nzg) tpr = fltarr(nlon,nlat,nzg,ntimes) gtpr = fltarr(ntimes,nzg) sigma_ped = fltarr(nlon,nlat,nzg,ntimes) gsigma_ped = fltarr(ntimes,nzg) hl_sigma_ped = fltarr(ntimes, nzg) hl_jheating = fltarr(ntimes, nzg) for n=0, ntimes-1 do begin for l=0, nzg-1 do begin for i = 0, nlon-1 do begin for j=0, nlat-1 do begin mind = min(abs(mzg(i,j,*,n)-GeometricHeight(l))) ind = where(abs(mzg(i,j,*,n)-GeometricHeight(l)) eq mind) ffp = ind(0) if (mzg(i,j,ffp,n) ge GeometricHeight(l)) then begin wa = (GeometricHeight(l)-mzg(i,j,ffp-1,n))/(mzg(i,j,ffp,n)-mzg(i,j,ffp-1,n)) wb = (mzg(i,j,ffp,n)-GeometricHeight(l))/(mzg(i,j,ffp,n)-mzg(i,j,ffp-1,n)) zg(i,j,l,n) = mzg(i,j,ffp-1,n)*wb+mzg(i,j,ffp,n)*wa den(i,j,l,n) = mden(i,j,ffp-1,n)*wb+mden(i,j,ffp,n)*wa eden(i,j,l,n) = meden(i,j,ffp-1,n)*wb+meden(i,j,ffp,n)*wa tn(i,j,l,n) = mtn(i,j,ffp-1,n)*wb+mtn(i,j,ffp,n)*wa tpr(i,j,l,n) = mtpr(i,j,ffp-1,n)*wb+mtpr(i,j,ffp,n)*wa sigma_ped(i,j,l,n) = msigma_ped(i,j,ffp-1,n)*wb+msigma_ped(i,j,ffp,n)*wa co2cool(i,j,l,n) = mco2cool(i,j,ffp-1,n)*wb+mco2cool(i,j,ffp,n)*wa nocool(i,j,l,n) = mnocool(i,j,ffp-1,n)*wb+mnocool(i,j,ffp,n)*wa sheating(i,j,l,n) = msheating(i,j,ffp-1,n)*den(i,j,l,n)*100.*wb+msheating(i,j,ffp,n)*den(i,j,l,n)*100.*wa ;mW/m3 sheating_pm(i,j,l,n) = msheating(i,j,ffp-1,n)*wb+msheating(i,j,ffp,n)*wa jheating(i,j,l,n) = mqjoule(i,j,ffp-1,n)*den(i,j,l,n)*100.*wb+mqjoule(i,j,ffp,n)*den(i,j,l,n)*100.*wa ;mW/m3 endif else begin wa = (mzg(i,j,ffp+1,n)-GeometricHeight(l))/(mzg(i,j,ffp+1,n)-mzg(i,j,ffp,n)) wb = (GeometricHeight(l)-mzg(i,j,ffp,n))/(mzg(i,j,ffp+1,n)-mzg(i,j,ffp,n)) zg(i,j,l,n) = mzg(i,j,ffp+1,n)*wb+mzg(i,j,ffp,n)*wa den(i,j,l,n) = mden(i,j,ffp+1,n)*wb+mden(i,j,ffp,n)*wa eden(i,j,l,n) = meden(i,j,ffp+1,n)*wb+meden(i,j,ffp,n)*wa tn(i,j,l,n) = mtn(i,j,ffp+1,n)*wb+mtn(i,j,ffp,n)*wa tpr(i,j,l,n) = mtpr(i,j,ffp+1,n)*wb+mtpr(i,j,ffp,n)*wa sigma_ped(i,j,l,n) = msigma_ped(i,j,ffp+1,n)*wb+msigma_ped(i,j,ffp,n)*wa co2cool(i,j,l,n) = mco2cool(i,j,ffp+1,n)*wb+mco2cool(i,j,ffp,n)*wa nocool(i,j,l,n) = mnocool(i,j,ffp+1,n)*wb+mnocool(i,j,ffp,n)*wa sheating(i,j,l,n) = msheating(i,j,ffp+1,n)*den(i,j,l,n)*100.*wb+msheating(i,j,ffp,n)*den(i,j,l,n)*100.*wa ;mW/m3 sheating_pm(i,j,l,n) = msheating(i,j,ffp+1,n)*wb+msheating(i,j,ffp,n)*wa jheating(i,j,l,n) = mqjoule(i,j,ffp+1,n)*den(i,j,l,n)*100.*wb+mqjoule(i,j,ffp,n)*den(i,j,l,n)*100.*wa ;mW/m3 endelse endfor ;j:lat endfor ;i:lon ; gzg(n,l) = mean(zg(*,*,l,n),/nan) ; gden(n,l) = mean(den(*,*,l,n),/nan) ; geden(n,l) = mean(eden(*,*,l,n),/nan) ; gtpr(n,l) = mean(tpr(*,*,l,n),/nan) ; gsigma_ped(n,l) = mean(sigma_ped(*,*,l,n),/nan) ; high_lat = where(abs(lat) ge 50) ; hl_sigma_ped(n,l) = mean(sigma_ped(*,high_lat,l,n),/nan) ; hl_jheating(n,l) = mean(jheating(*,high_lat,l,n),/nan) ;mW/m3 ; gtn(n,l) = mean(tn(*,*,l,n),/nan) ; gden(n,l) = mean(den(*,*,l,n),/nan) ; gco2cool(n,l) = mean(co2cool(*,*,l,n)*1.e-4,/nan) ;mW/g ; gnocool(n,l) = mean(nocool(*,*,l,n)*1.e-4,/nan) ;mW/g ; gsheating(n,l) = mean(sheating(*,*,l,n),/nan) ;mW/m3 ; gmsheating(n,l) = mean(sheating_pm(*,*,l,n),/nan)*1.e-4 ;mW/g ; gjheating(n,l) = mean(jheating(*,*,l,n),/nan) ;mW/m3 endfor ;l endfor ;ntimes ;file 22222222222222222222222222222222222-------------------------------------------------------------- ;files2 = file_search('$LINUX/SolarFlare/nc2003301/sfism_0_195nm_2003301.nc') files2 = file_search('/hao/aim3/huangys/tiegcm1.94_data/SolarFlare/nc2003301/sfism_0_14nm_2003301.nc') fn = '2003301_0_14nm' nfiles2 = n_elements(files2) print, 'number of files = ',nfiles2 mtime2 = intarr(3,288*nfiles2) mut2 = fltarr(288*nfiles2) mhpower2 = fltarr(288*nfiles2) msigma_ped2 = fltarr(72,36,29,288*nfiles2);S/m meden2 = fltarr(72,36,29,288*nfiles2);ilev mden2 = fltarr(72,36,29,288*nfiles2);ilev mqjoule2 = fltarr(72,36,29,288*nfiles2) ; (lon.lat,lev,time) mqjoule_integ2 = fltarr(72,36,288*nfiles2); (lon,lat,time) msheating2 = fltarr(72,36,29,288*nfiles2) mtn2 = fltarr(72,36,29,288*nfiles2) mzg2 = fltarr(72,36,29,288*nfiles2);ilev mz2 = fltarr(72,36,29,288*nfiles2);ilev mtec2 = fltarr(72,36,288*nfiles2) mqop2pa2 = fltarr(72,36,29,288*nfiles2) ; production rate (cm-3s-1) mqop2da2 = fltarr(72,36,29,288*nfiles2) ; (lon,lat,lev,time) mqo2p2 = fltarr(72,36,29,288*nfiles2) ; (lon,lat,lev,time) mqn2p2 = fltarr(72,36,29,288*nfiles2) ; (lon,lat,lev,time) mqnp2 = fltarr(72,36,29,288*nfiles2) ; (lon,lat,lev,time) mqopa2 = fltarr(72,36,29,288*nfiles2) ; (lon,lat,lev,time) mqnop2 = fltarr(72,36,29,288*nfiles2) ; (lon,lat,lev,time) mco2cool2 = fltarr(72,36,29,288*nfiles) mnocool2 = fltarr(72,36,29,288*nfiles) for i = 0L,nfiles2-1 do begin print, 'i=',i filein2 = files2(i) ;; open and read netCDF file print, 'Reading ',filein2 id=ncdf_open(filein2) ncdf_varget,id,'mtime',var mtime2(*,i*288:(i*288+287)) = var ncdf_varget,id,ncdf_varid(id,'ut'),vut ;universal time (hours) mut2(i*288:(i*288+287)) = vut ncdf_varget,id,ncdf_varid(id,'CO2_COOL'),vco2 ;erg/g/s mco2cool2(*,*,*,i*288:(i*288+287)) = vco2 ncdf_varget,id,ncdf_varid(id,'NO_COOL'),vno ;erg/g/s mnocool2(*,*,*,i*288:(i*288+287)) = vno ncdf_varget,id,ncdf_varid(id,'hpower'),vhpower ;hemispheric power, GW mhpower2(i*288:(i*288+287)) = vhpower ncdf_varget,id,ncdf_varid(id,'SIGMA_PED'),vped ;pederson conductivity msigma_ped2(*,*,*,i*288:(i*288+287)) = vped ncdf_varget,id,ncdf_varid(id,'NE'),veden ;total density, g/cm3 meden2(*,*,*,i*288:(i*288+287)) = veden ncdf_varget,id,ncdf_varid(id,'DEN'),vden ;total density, g/cm3 mden2(*,*,*,i*288:(i*288+287)) = vden ncdf_varget,id,ncdf_varid(id,'QJOULE'),vqj mqjoule2(*,*,*,i*288:(i*288+287)) = vqj ncdf_varget,id,ncdf_varid(id,'QJOULE_INTEG'),vqj_int ;erg/cm2/s = mW/m2 mqjoule_integ2(*,*,i*288:(i*288+287)) = vqj_int ncdf_varget,id,ncdf_varid(id,'SHEATING'),vheat ;solar heating,erg/g/s msheating2(*,*,*,i*288:(i*288+287)) = vheat ncdf_varget,id,ncdf_varid(id,'TN'),vtn mtn2(*,*,*,i*288:(i*288+287)) = vtn ncdf_varget,id,ncdf_varid(id,'ZG'),vzg mzg2(*,*,*,i*288:(i*288+287)) = vzg ncdf_varget,id,ncdf_varid(id,'Z'),vz mz2(*,*,*,i*288:(i*288+287)) = vz ncdf_varget,id,ncdf_varid(id,'TEC'),vtec mtec2(*,*,i*288:(i*288+287)) = vtec ;1/cm2 ncdf_varget,id,ncdf_varid(id,'QOP2Pa'),vqop2pa mqop2pa2(*,*,*,i*288:(i*288+287)) = vqop2pa ncdf_varget,id,ncdf_varid(id,'QOP2Da'),vqop2da mqop2da2(*,*,*,i*288:(i*288+287)) = vqop2da ncdf_varget,id,ncdf_varid(id,'QO2P'),vqo2p mqo2p2(*,*,*,i*288:(i*288+287)) = vqo2p ncdf_varget,id,ncdf_varid(id,'QN2P'),vqn2p mqn2p2(*,*,*,i*288:(i*288+287)) = vqn2p ncdf_varget,id,ncdf_varid(id,'QNP'),vqnp mqnp2(*,*,*,i*288:(i*288+287)) = vqnp ncdf_varget,id,ncdf_varid(id,'QOPa'),vqopa mqopa2(*,*,*,i*288:(i*288+287)) = vqopa ncdf_varget,id,ncdf_varid(id,'QNOP'),vqnop mqnop2(*,*,*,i*288:(i*288+287)) = vqnop print,'Reading Finished~~~~~~~~~~~~~~~' endfor mtime2 = transpose(mtime2) ntimes2 = n_elements(mtime2(*,0)) ;# of hours smtime2 = fltarr(ntimes2) smtime2(*) = mtime2(*,0)+mtime2(*,1)/24.+mtime2(*,2)/1440. fp = where(mco2cool2 ge missing) mco2cool2(fp) = !Values.F_NAN ;NaN fp = where(mnocool2 ge missing) mnocool2(fp) = !Values.F_NAN ;NaN fp = where(mhpower2 ge missing) mhpower2(fp) = !Values.F_NAN ;NaN fp = where(msigma_ped2 ge missing) msigma_ped2(fp) = !Values.F_NAN ;NaN fp = where(meden2 ge missing) meden2(fp) = !Values.F_NAN ;NaN fp = where(mden2 ge missing) mden2(fp) = !Values.F_NAN ;NaN fp = where(mqjoule2 ge missing) mqjoule2(fp) = !Values.F_NAN ;NaN fp = where(mqjoule_integ2 ge missing) mqjoule_integ2(fp) = !Values.F_NAN ;NaN fp = where(msheating2 ge missing) msheating2(fp) = !Values.F_NAN ;NaN fp = where(mtn2 ge missing) mtn2(fp) = !Values.F_NAN ;NaN fp = where(mzg2 ge missing) mzg2(fp) = !Values.F_NAN ;NaN fp = where(mz2 ge missing) mz2(fp) = !Values.F_NAN ;NaN fp = where(mtec2 ge missing) mtec2(fp) = !Values.F_NAN ;NaN fp = where(mqop2pa2 ge missing) mqop2pa2(fp) = !Values.F_NAN ;NaN fp = where(mqop2da2 ge missing) mqop2da2(fp) = !Values.F_NAN ;NaN fp = where(mqo2p2 ge missing) mqo2p2(fp) = !Values.F_NAN ;NaN fp = where(mqn2p2 ge missing) mqn2p2(fp) = !Values.F_NAN ;NaN fp = where(mqnp2 ge missing) mqnp2(fp) = !Values.F_NAN ;NaN fp = where(mqopa2 ge missing) mqopa2(fp) = !Values.F_NAN fp = where(mqnop2 ge missing) mqnop2(fp) = !Values.F_NAN ;NaN mtpr2 = mqop2pa2+mqop2da2+mqo2p2+mqn2p2+mqnp2+mqopa2+mqnop2 ;total production rate ;-------------------------Finish loading data tn2 = fltarr(nlon,nlat,nzg,ntimes2) gtn2 = fltarr(ntimes2,nzg) zg2 = fltarr(nlon,nlat,nzg,ntimes2) gzg2 = fltarr(ntimes2,nzg) den2 = fltarr(nlon,nlat,nzg,ntimes2) gden2 = fltarr(ntimes2,nzg) eden2 = fltarr(nlon,nlat,nzg,ntimes2) geden2 = fltarr(ntimes2,nzg) co2cool2 = fltarr(nlon,nlat,nzg,ntimes2) gco2cool2 = fltarr(ntimes2,nzg) nocool2 = fltarr(nlon,nlat,nzg,ntimes2) gnocool2 = fltarr(ntimes2,nzg) sheating2=fltarr(nlon,nlat,nzg,ntimes2) sheating2_pm=fltarr(nlon,nlat,nzg,ntimes2) ; solar heating per mass gsheating2=fltarr(ntimes2,nzg) gmsheating2=fltarr(ntimes2,nzg) jheating2=fltarr(nlon,nlat,nzg,ntimes2) gjheating2=fltarr(ntimes2,nzg) tpr2 = fltarr(nlon,nlat,nzg,ntimes2) gtpr2 = fltarr(ntimes2,nzg) sigma_ped2 = fltarr(nlon,nlat,nzg,ntimes2) gsigma_ped2 = fltarr(ntimes2,nzg) hl_sigma_ped2 = fltarr(ntimes2, nzg) hl_jheating2 = fltarr(ntimes2, nzg) for n=0, ntimes2-1 do begin for l=0, nzg-1 do begin for i = 0, nlon-1 do begin for j=0, nlat-1 do begin mind = min(abs(mzg2(i,j,*,n)-GeometricHeight(l))) ind = where(abs(mzg2(i,j,*,n)-GeometricHeight(l)) eq mind) ffp = ind(0) if (mzg2(i,j,ffp,n) ge GeometricHeight(l)) then begin wa = (GeometricHeight(l)-mzg2(i,j,ffp-1,n))/(mzg2(i,j,ffp,n)-mzg2(i,j,ffp-1,n)) wb = (mzg2(i,j,ffp,n)-GeometricHeight(l))/(mzg2(i,j,ffp,n)-mzg2(i,j,ffp-1,n)) zg2(i,j,l,n) = mzg2(i,j,ffp-1,n)*wb+mzg2(i,j,ffp,n)*wa den2(i,j,l,n) = mden2(i,j,ffp-1,n)*wb+mden2(i,j,ffp,n)*wa eden2(i,j,l,n) = meden2(i,j,ffp-1,n)*wb+meden2(i,j,ffp,n)*wa tn2(i,j,l,n) = mtn2(i,j,ffp-1,n)*wb+mtn2(i,j,ffp,n)*wa tpr2(i,j,l,n) = mtpr2(i,j,ffp-1,n)*wb+mtpr2(i,j,ffp,n)*wa sigma_ped2(i,j,l,n) = msigma_ped2(i,j,ffp-1,n)*wb+msigma_ped2(i,j,ffp,n)*wa co2cool2(i,j,l,n) = mco2cool2(i,j,ffp-1,n)*wb+mco2cool2(i,j,ffp,n)*wa nocool2(i,j,l,n) = mnocool2(i,j,ffp-1,n)*wb+mnocool2(i,j,ffp,n)*wa sheating2(i,j,l,n) = msheating2(i,j,ffp-1,n)*den2(i,j,l,n)*100.*wb+msheating2(i,j,ffp,n)*den2(i,j,l,n)*100.*wa ;mW/m3 sheating2_pm(i,j,l,n) = msheating2(i,j,ffp-1,n)*wb+msheating2(i,j,ffp,n)*wa jheating2(i,j,l,n) = mqjoule2(i,j,ffp-1,n)*den2(i,j,l,n)*100.*wb+mqjoule2(i,j,ffp,n)*den2(i,j,l,n)*100.*wa ;mW/m3 endif else begin wa = (mzg2(i,j,ffp+1,n)-GeometricHeight(l))/(mzg2(i,j,ffp+1,n)-mzg2(i,j,ffp,n)) wb = (GeometricHeight(l)-mzg2(i,j,ffp,n))/(mzg2(i,j,ffp+1,n)-mzg2(i,j,ffp,n)) zg2(i,j,l,n) = mzg2(i,j,ffp+1,n)*wb+mzg2(i,j,ffp,n)*wa den2(i,j,l,n) = mden2(i,j,ffp+1,n)*wb+mden2(i,j,ffp,n)*wa eden2(i,j,l,n) = meden2(i,j,ffp+1,n)*wb+meden2(i,j,ffp,n)*wa tn2(i,j,l,n) = mtn2(i,j,ffp+1,n)*wb+mtn2(i,j,ffp,n)*wa tpr2(i,j,l,n) = mtpr2(i,j,ffp+1,n)*wb+mtpr2(i,j,ffp,n)*wa sigma_ped2(i,j,l,n) = msigma_ped2(i,j,ffp+1,n)*wb+msigma_ped2(i,j,ffp,n)*wa co2cool2(i,j,l,n) = mco2cool2(i,j,ffp+1,n)*wb+mco2cool2(i,j,ffp,n)*wa nocool2(i,j,l,n) = mnocool2(i,j,ffp+1,n)*wb+mnocool2(i,j,ffp,n)*wa sheating2(i,j,l,n) = msheating2(i,j,ffp+1,n)*den2(i,j,l,n)*100.*wb+msheating2(i,j,ffp,n)*den2(i,j,l,n)*100.*wa ;mW/m3 sheating2_pm(i,j,l,n) = msheating2(i,j,ffp+1,n)*wb+msheating2(i,j,ffp,n)*wa jheating2(i,j,l,n) = mqjoule2(i,j,ffp+1,n)*den2(i,j,l,n)*100.*wb+mqjoule2(i,j,ffp,n)*den2(i,j,l,n)*100.*wa ;mW/m3 endelse endfor ;j:lat endfor ;i:lon ; gzg2(n,l) = mean(zg2(*,*,l,n),/nan) ; gden2(n,l) = mean(den2(*,*,l,n),/nan) ; geden2(n,l) = mean(eden2(*,*,l,n),/nan) ; gtpr2(n,l) = mean(tpr2(*,*,l,n),/nan) ; gsigma_ped2(n,l) = mean(sigma_ped2(*,*,l,n),/nan) ; high_lat = where(abs(lat) ge 50) ; hl_sigma_ped2(n,l) = mean(sigma_ped2(*,high_lat,l,n),/nan) ; hl_jheating2(n,l) = mean(jheating2(*,high_lat,l,n),/nan) ;mW/m3 ; gtn2(n,l) = mean(tn2(*,*,l,n),/nan) ; gden2(n,l) = mean(den2(*,*,l,n),/nan) ; gco2cool2(n,l) = mean(co2cool2(*,*,l,n)*1.e-4,/nan);mW/g ; gnocool2(n,l) = mean(nocool2(*,*,l,n)*1.e-4,/nan);mW/g ; gsheating2(n,l) = mean(sheating2(*,*,l,n),/nan) ;mW/m3 ; gmsheating2(n,l) = mean(sheating2_pm(*,*,l,n),/nan)*1.e-4 ;mW/g ; gjheating2(n,l) = mean(jheating2(*,*,l,n),/nan) ;mW/m3 endfor ;l endfor ;ntimes ;calculate the percentage increase of No cooling between preflare and XUV perc_nocool = (nocool2-nocool)/nocool*100. ;(lon,lat,1,time) ;print, size(perc_nocool) diff_nocool = nocool2-nocool sample_time = (180./60.+11.)/24.+301. ;flare starts at 11UT min_time = min(abs(smtime-sample_time), fp_time) result = reform(perc_nocool(*,*,0,fp_time)) result2 = reform(diff_nocool(*,*,0,fp_time)) print, smtime(fp_time) lon = lon+180. set_plot,'ps' device,filename=fn+'_nocool_180min.eps', /color conmake,result,lon,lat,256,/cbar,xran=[0, 360],yran=[-90,90],title = 'NOcooling %increase '+fn, $ leve=findgen(256)/255*50.1,minlev = 0, maxlev = 50.5 conmake,result2,lon,lat,256,/cbar,xran=[0, 360],yran=[-90,90],title = 'NOcooling increase(erg/g/s) '+fn, $ leve=findgen(256)/255*5001.,minlev = 0., maxlev = 5005. device,/close set_plot,'x' stop ;difference(n,l)-------------------------------------------------------- dgco2cool = gco2cool2-gco2cool dgnocool = gnocool2-gnocool dgjheating=gjheating2-gjheating dgsheating=gsheating2-gsheating dgmsheating=gmsheating2-gmsheating dgeden=geden2-geden dgden=gden2-gden dgtn=gtn2-gtn dgped=gsigma_ped2-gsigma_ped dgtpr = gtpr2-gtpr dhljheating = hl_jheating2 - hl_jheating dhlped = hl_sigma_ped2 - hl_sigma_ped ;difference(lon,lat,lev,time)------------------------------------------- dco2cool = mco2cool2-mco2cool dnocool = mnocool2-mnocool ;set_plot,'ps' ;device,decomposed=0,/color ;device,filename='dcooling.ps' loadct, 33 number_levels = 33 levels_contour = dblarr(number_levels) colors_contour = dblarr(number_levels) color_table_offset = 2.0 ;tind = where(abs(lat) le 50.0) maxData = max(dco2cool(*,*,20,180)) minData = min(dco2cool(*,*,20,180)) ;!p.multi = [0,2,1] for i = 0, number_levels-1 do begin levels_contour(i) = (maxData-minData)*i/(number_levels-1.0)+minData colors_contour(i) = (255.0-color_table_offset)*i/(number_levels-1.0)+color_table_offset endfor contour, dco2cool(*,*,20,180),lon,lat, $ charsize=1, charthick=2, levels=levels_contour,c_color=colors_contour, $ xrange=[-180,180],/fill,position=[0.2,0.4,0.9,0.9], $ xstyle=1,ystyle=1, $ xticks=7,xtickv=[-180,-120,-60,0,60,120,180],xtickname=['0','60','120','180','240','300','360'];,/overplot contour, [[levels_contour],[levels_contour]],levels_contour,indgen(2),charsize=1,charthick=2, $ levels=levels_contour,c_color=colors_contour,/fill,position=[0.2,0.25,0.9,0.28],/noerase,ystyle=4, $ xr=[minData,maxData],xstyle=1;,XTICKFORMAT="(A1)" maxData = max(dnocool(*,*,20,180)) minData = min(dnocool(*,*,20,180)) for i = 0, number_levels-1 do begin levels_contour(i) = (maxData-minData)*i/(number_levels-1.0)+minData colors_contour(i) = (255.0-color_table_offset)*i/(number_levels-1.0)+color_table_offset endfor contour, dnocool(*,*,20,180),lon,lat, $ charsize=1, charthick=2, levels=levels_contour,c_color=colors_contour, $ xrange=[-180,180],/fill,position=[0.2,0.4,0.9,0.9], $ xstyle=1,ystyle=1, $ xticks=7,xtickv=[-180,-120,-60,0,60,120,180],xtickname=['0','60','120','180','240','300','360'];,/overplot contour, [[levels_contour],[levels_contour]],levels_contour,indgen(2),charsize=1,charthick=2, $ levels=levels_contour,c_color=colors_contour,/fill,position=[0.2,0.25,0.9,0.28],/noerase,ystyle=4, $ xr=[minData,maxData],xstyle=1;,XTICKFORMAT="(A1)" ;dd = transpose(reform(den2(0,*,20,*)-den(0,*,20,*))) ;maxData = max(dd) ;minData = min(dd) ;for i = 0, number_levels-1 do begin ; levels_contour(i) = (maxData-minData)*i/(number_levels-1.0)+minData ; colors_contour(i) = (255.0-color_table_offset)*i/(number_levels-1.0)+color_table_offset ;endfor ;contour,dd,smtime,lat, $ ;charsize=1, charthick=2, levels=levels_contour,c_color=colors_contour,/fill,position=[0.2,0.4,0.9,0.9] ;contour, [[levels_contour],[levels_contour]],levels_contour,indgen(2),charsize=1,charthick=2, $ ;levels=levels_contour,c_color=colors_contour,/fill,position=[0.2,0.25,0.9,0.28],/noerase,ystyle=4, $ ;xr=[minData,maxData],xstyle=1;,XTICKFORMAT="(A1)" ;xyouts, -5,-3,'MIN = '+strmid(strtrim(min(up(*,tind)),2),0,5) ;xyouts, 0,-3,'MAX = '+strmid(strtrim(max(up(*,tind)),2),0,4) ;device,/close ;set_plot,'x' ;stop TI_sh = total(dgsheating,1,/nan) TI_jh = total(dgjheating,1,/nan) hlTI_jh = total(dhljheating,1,/nan) shp = dgsheating/gsheating*100. mshp = dgmsheating/gmsheating*100. jhp = dgjheating/gjheating*100. hl_jhp = dhljheating/hl_jheating*100. TI_shp = total(shp, 1, /nan) TI_jhp = total(jhp, 1, /nan) ;;maximum perturbation in electron density: dgeden(288, 29)-------------------------------------- ;max_eden = max(dgeden,i,/nan) ;max_eden_smtime = i mod 288 ;max_eden_lev = i/288 ;max_ped = max(dgped,i,/nan) ;max_ped_smtime = i mod 288 ;max_jh = max(dgjheating,i,/nan) ;max_jh_smtime = i mod 288 ;max_sh = max(dgsheating,i,/nan) ;max_sh_smtime = i mod 288 ;print, smtime(max_ped_smtime),smtime(max_ped_smtime),smtime(max_jh_smtime),smtime(max_sh_smtime) ;!p.multi = [0,1,2] ;plot, dgeden(max_eden_smtime,*),lev ;plot, dgped(max_ped_smtime,*),lev ;!p.multi =0 den_pct = dgden/gden*100 max_den_pct = fltarr(nzg) max_tn = fltarr(nzg) ptime_den = fltarr(nzg) ptime_tn = fltarr(nzg) for l=3, nzg-1 do begin; higher than 145km max_den_pct(l) = max(den_pct(*,l), location,/nan) ptime_den(l) = smtime(location) max_tn(l) = max(dgtn(*,l), location,/nan) ptime_tn(l) = smtime(location) ;;ind = array_indices(dgden/gden*100,location) ;;print, smtime(ind(0)),GeometricHeight(ind(1)) endfor ;set_plot,'ps' ;device,filename='Peaktime'+fn+'.eps',/portrait ;plot, ptime_tn, GeometricHeight/1.e5 ;xyouts, 301.2, 400, 'Tn' ;plot, ptime_den, GeometricHeight/1.e5 ;xyouts, 301.2, 400, 'Den' ;device,/close ;set_plot,'x' timeind = 134 ;180 for smtime = 301.625 and 134 for smtime at 301.465 (at peak) ;openw,33,'hl_dgeden_zg'+fn+'_3016.dat',width = 300 openw,33,'dgmsheating_zg_peak'+fn+'.dat',width = 300 ;printf,33,'zg(km) dgeden(cm-3) dgped dgjheating dgmsheating ;dgtpr(cm-3s-1) dgden dgtn TI_jh TI_sh TI_jhp TI_shp hlTI_jh high_lat_jh high_lat_ped' for n = 0, nzg-1 do begin printf,33,GeometricHeight(n),dgeden(timeind,n),dgped(timeind,n),dgjheating(timeind,n),dgmsheating(timeind,n),dgtpr(timeind, n),dgden(timeind,n),dgtn(timeind,n),TI_jh(n),TI_sh(n),TI_jhp(n), TI_shp(n), hlTI_jh(n), dhljheating(timeind,n), dhlped(timeind,n) endfor close,33 timeind = 180 ;for smtime = 301.625 ;openw,33,'hl_dgeden_zg'+fn+'_3016.dat',width = 300 openw,35,'dgeden_zg_3pm'+fn+'.dat',width = 300 ;printf,35,'zg(km) dgeden(cm-3) dgped dgjheating dgsheating ;dgtpr(cm-3s-1) dgden dgtn TI_jh TI_sh TI_jhp TI_shp hlTI_jh high_lat_jh high_lat_ped' for n = 0, nzg-1 do begin printf,35,GeometricHeight(n),dgeden(timeind,n),dgped(timeind,n),dgjheating(timeind,n),dgmsheating(timeind,n),dgtpr(timeind, n),dgden(timeind,n),dgtn(timeind,n),TI_jh(n),TI_sh(n),TI_jhp(n), TI_shp(n), hlTI_jh(n), dhljheating(timeind,n), dhlped(timeind,n) endfor close,35 ;stop l_400 =20 openw,13,'dgmsheating_400_time'+fn+'.dat',width=300 printf,13,'smtime(day) dgmsheating_400(mW/g) dgjheating_400(mW/m3) dgden_400(g/cm3) dgtn_400(K) dgped_400(S/m) dgtpr_400(cm-3s-1) dgeden_400(cm-3) dhljheating_400(mW/m3) dhlped_400(mW/m3)' for n = 0, n_elements(smtime)-1 do begin printf,13,smtime(n), dgmsheating(n,l_400),dgjheating(n,l_400),dgden(n,l_400),dgtn(n,l_400),dgped(n,l_400),dgtpr(n,l_400),dgeden(n,l_400),dhljheating(n,l_400),dhlped(n,l_400) endfor close,13 ;stop ;!p.multi = [0,1,2] set_plot,'ps' ;device,filename=fn+'_zg_3.eps', /color;/ENCAPSUL, BITS_PER_PIXEL=8;, /COLOR device,filename=fn+'_zg_pm_5.eps', /color;/ENCAPSUL, BITS_PER_PIXEL=8;, /COLOR xrange = [301,302] conmake,dgmsheating,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'GSH(mW/g)'+fn;, overp=dgsheating ;plot, smtime, dgsheating(*,l_400),xran=xrange,ytitle='GSH at 400 km(mW/m3)',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 ;conmake,dhljheating,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'HLJH(mw/m3)'+fn;, overp=dgjheating ;plot, smtime, dhljheating(*,l_400),xran=xrange,ytitle='HLJH at 400 km(mW/m3)',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 conmake,dgden,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title= 'Gdensity (g/cm3)'+fn;,minlev=-2.0e-14,maxlev =2.0e-14;, overp=dgden ;plot, smtime, dgden(*,l_400),xran=xrange,ytitle='Gdensity at 400 km(g/cm3)',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 conmake,dgtn,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'Gtemperature(K)'+fn,leve=findgen(256)/255*30., oxx=ptime_tn(3:29), oyy=GeometricHeight(3:29)/1.e5;,minlev=0.0 ;loadct,39 ;oplot, ptime_tn, GeometricHeight/1.e5,thick = 2,color = 255;!white ;plot, smtime, dgtn(*,l_400),xran=xrange,ytitle='Gtemperature at 400 km(K)',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 ;conmake,dhlped, smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'HLPederson(S/m)'+fn;, overp=dgtn ;plot, smtime, dhlped(*,l_400),xran=xrange,ytitle='HLPederson 400 km(S/m)',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 ;conmake,dgtpr, smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'GProduction rate(cm-3s-1)'+fn ;plot, smtime, dgtpr(*,l_400),xran=xrange,ytitle='GProduction rate at 400 km(cm-3s-1)',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 conmake,dgeden, smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'GNe(cm-3)'+fn,leve=findgen(256)/255*1.e5;,minlev=0.0 ;plot, smtime, dgeden(*,l_400),xran=xrange,ytitle='GNe at 400 km(cm-3)',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 conmake,dgco2cool,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'GCO2cooling(mW/g)'+fn conmake,dgnocool,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'GNOcoolingy(mW/g)'+fn device,/close set_plot,'x' ;!p.multi = [0,1,2] set_plot,'ps' ;device,filename=fn+'_zg_4.eps', /color;/ENCAPSUL, BITS_PER_PIXEL=8;, /COLOR device,filename=fn+'_zg_pm_6.eps', /color;/ENCAPSUL, BITS_PER_PIXEL=8;, /COLOR xrange = [301,302] conmake,mshp,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100, 500],title = '0-175 nm',leve=findgen(256)/255*100.;,minlev=0.0 ;plot, smtime, dgsheating(*,l_400)/gsheating(*,l_400)*100.,xran=xrange,ytitle='GSH %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 ;conmake,hl_jhp,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100, 500],title = 'HLJH %increase'+fn;,maxlev =30.;, overp=dgjheating ;plot, smtime, dhljheating(*,l_400)/hl_jheating(*,l_400)*100.,xran=xrange,ytitle='HLJH %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 conmake,dgden/gden*100.,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100, 500],title= 'Gdensity %increase'+fn,leve=findgen(256)/255*10.,oxx=ptime_den(3:29),oyy=GeometricHeight(3:29)/1.e5;,minlev=0.0 ;loadct, 39 ;oplot, ptime_den, GeometricHeight/1.e5,thick = 2,color = 255;!white ;plot, smtime, dgden(*,l_400)/gden(*,l_400)*100.,xran=xrange,ytitle='Gdensity %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 conmake,dgtn/gtn*100.,smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100, 500],title = 'Gtemperature %increase'+fn;, overp=dgtn ;plot, smtime, dgtn(*,l_400)/gtn(*,l_400)*100.,xran=xrange,ytitle='Gtemperature %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 ;conmake,dhlped/hl_sigma_ped*100., smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100, 500],title = 'HLPederson %increase'+fn;,maxlev=30.;, overp=dgtn ;plot, smtime, dhlped(*,l_400)/hl_sigma_ped(*,l_400)*100.,xran=xrange,ytitle='HLPederson %increase 400 km',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 ;conmake,dgtpr/gtpr*100., smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100, 500],title = 'GProduction rate %increase'+fn ;plot, smtime, dgtpr(*,l_400)/gtpr(*,l_400)*100.,xran=xrange,ytitle='GProduction rate %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 conmake,dgeden/geden*100., smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100, 500],title = 'GNe %increase'+fn ;plot, smtime, dgeden(*,l_400)/geden(*,l_400)*100.,xran=xrange,ytitle='GNe %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]] ;oplot, [250,311],[0,0],linestyle=2 conmake,dgco2cool/gco2cool*100, smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'GCO2cooling %increase'+fn conmake,dgnocool/gnocool*100, smtime,GeometricHeight/1.e5,256,/cbar,xlabel=5,ylabel=1,xran=xrange,yran=[100,500],title = 'GNOcooling %increase'+fn device,/close set_plot,'x' stop ;total integration------------------------------------ ;xrange = [301.,302.] !p.multi=[0,1,2] set_plot,'ps' device,filename='Energyinputs_'+fn+'_3.eps',/portrait ;plot, smtime2,int_sh2,ytitle='solar flux (GW)',title = fn ;oplot, smtime,int_sh,linestyle=2 plot, smtime, dsh,ytitle='dSH (GW)',title=fn,xran=xrange xyouts, 301.7, max(dsh,/nan)/2.,string(total_sh)+'J' plot, smtime, djh,ytitle='dJH (GW)',xran=xrange xyouts, 301.7, max(djh,/nan)/2., string(total_jh)+'J' device,/close set_plot,'x' stop ;polar contour---------------------------------------------------------- !p.multi= 0 lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180. mqjoule_integ_t = fltarr(nlon+1,nlat,ntimes) & mqjoule_integ_t(0:71,*,*)= mqjoule_integ & mqjoule_integ_t(72,*,*)=mqjoule_integ(0,*,*) set_plot,'ps' device,filename='polar.eps',/color nl =256 loadct,34, ncolors = nl ;LOADCT, 33,ncolors = nl DEVICE, DECOMPOSED = 0 map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[50.,-180,90.,180] contour,mqjoule_integ_t(*,18:35,150),lon2,lat(18:35),$ nlevels=nl,C_COLORS = INDGEN(nl),xstyle=1,ystyle=1, $ /overplot,/fill,/isotropic,POSITION=[ 0.14, 0.08, 0.82, 0.92 ] ;,position=position,zrange=[0.,2.6],zstyle= 1,c_colors = c_colors map_grid,/label,latlab = 175,color=130 map_continents,color=130 ;color bar mxx = MAX( mqjoule_integ_t, MIN = mnn ,/nan) * 0.98 delta_lev = ( mxx - mnn ) / FLOAT( nl ) level = FINDGEN( nl ) * delta_lev + mnn ;colorbar array consists of levels previously determined cbar_arr = FLTARR( 2, N_ELEMENTS( level ) ) cbar_arr[0,*] = level cbar_arr[1,*] = level CONTOUR, cbar_arr, [1,2], level, $ NLEVELS = nl, C_COLORS = INDGEN( nl ), XSTYLE = 5, YSTYLE = 5, $ /ZSTYLE, /FILL, BACKGROUND = nl - 1, COLOR = 0, $ POSITION = [ 0.82 + 0.03, 0.03, $ 0.82 + 0.06, 0.945],$ LEVELS = level, /NOERASE, XTICKS = 1, XTICKNAME = [' ', ' ', ' '] AXIS, /YAXIS, charsize=1., charthick=2, /DATA, COLOR = 0, $ /ZSTYLE, /YSTYLE device,/close set_plot,'x' lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180. height_sh_t = fltarr(nlon+1,nlat,ntimes) & height_sh_t(0:71,*,*)= height_sh & height_sh_t(72,*,*)=height_sh(0,*,*) set_plot,'ps' device,filename='polar_sh.eps',/color nl =256 LOADCT,34,ncolors = nl DEVICE, DECOMPOSED = 0 map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[50.,-180,90.,180] contour,height_sh_t(*,18:35,0),lon2,lat(18:35),$ nlevels=nl,C_COLORS = INDGEN(nl),xstyle=1,ystyle=1, $ /overplot,/fill,/isotropic,POSITION=[ 0.14, 0.08, 0.82, 0.92 ] ;,position=position,zrange=[0.,2.6],zstyle= 1,c_colors = c_colors map_grid,/label,latlab = 175,color=130 map_continents,color=130 ;color bar mxx = MAX(height_sh, MIN = mnn ,/nan) * 0.98 delta_lev = ( mxx - mnn ) / FLOAT( nl ) level = FINDGEN( nl ) * delta_lev + mnn ;colorbar array consists of levels previously determined cbar_arr = FLTARR( 2, N_ELEMENTS( level ) ) cbar_arr[0,*] = level cbar_arr[1,*] = level CONTOUR, cbar_arr, [1,2], level, $ NLEVELS = nl, C_COLORS = INDGEN( nl ), XSTYLE = 5, YSTYLE = 5, $ /ZSTYLE, /FILL, BACKGROUND = nl - 1, COLOR = 0, $ POSITION = [ 0.82 + 0.03, 0.03, $ 0.82 + 0.06, 0.945],$ LEVELS = level, /NOERASE, XTICKS = 1, XTICKNAME = [' ', ' ', ' '] AXIS, /YAXIS, charsize=1., charthick=2, /DATA, COLOR = 0, $ /ZSTYLE, /YSTYLE device,/close set_plot,'x' lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180. tec_t = fltarr(nlon+1,nlat,ntimes) & tec_t(0:71,*,*)= mtec & tec_t(72,*,*)=mtec(0,*,*) tec_t2 = fltarr(nlon+1,nlat,ntimes2) & tec_t2(0:71,*,*)= mtec2 & tec_t2(72,*,*)=mtec2(0,*,*) dtec_t = tec_t-tec_t2 set_plot,'ps' device,filename='polar_tec.eps',/color nl =256 LOADCT, 34,ncolors = nl DEVICE, DECOMPOSED = 0 map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[50.,-180,90.,180] contour,dtec_t(*,*,150),lon2,lat,$ nlevels=nl,C_COLORS = INDGEN(nl),xstyle=1,ystyle=1, $ /overplot,/fill,/isotropic,POSITION=[ 0.14, 0.08, 0.82, 0.92 ] ;,position=position,zrange=[0.,2.6],zstyle= 1,c_colors = c_colors map_grid,/label,latlab = 175,color=130 map_continents,color=130 ;color bar mxx = MAX(dtec_t(*,*,150), MIN = mnn ,/nan) * 0.98 delta_lev = ( mxx - mnn ) / FLOAT( nl ) level = FINDGEN( nl ) * delta_lev + mnn ;colorbar array consists of levels previously determined cbar_arr = FLTARR( 2, N_ELEMENTS( level ) ) cbar_arr[0,*] = level cbar_arr[1,*] = level CONTOUR, cbar_arr, [1,2], level, $ NLEVELS = nl, C_COLORS = INDGEN( nl ), XSTYLE = 5, YSTYLE = 5, $ /ZSTYLE, /FILL, BACKGROUND = nl - 1, COLOR = 0, $ POSITION = [ 0.82 + 0.03, 0.03, $ 0.82 + 0.06, 0.945],$ LEVELS = level, /NOERASE, XTICKS = 1, XTICKNAME = [' ', ' ', ' '] AXIS, /YAXIS, charsize=1., charthick=2, /DATA, COLOR = 0, $ /ZSTYLE, /YSTYLE device,/close set_plot,'x' lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180. ped2 = fltarr(nlon+1,nlat) & ped2(0:71,*)= msigma_ped2(*,*,10,134) & ped2(72,*)=msigma_ped2(0,*,10,134) ped = fltarr(nlon+1,nlat) & ped(0:71,*)= msigma_ped(*,*,10,134) & ped(72,*)=msigma_ped(0,*,10,134) set_plot,'ps' device,filename='polar_pedersen.eps',/color nl =256 LOADCT, 34,ncolors = nl DEVICE, DECOMPOSED = 0 map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[30.,-180,90.,180] contour,ped2-ped,lon2,lat,$ nlevels=nl,C_COLORS = INDGEN(nl),xstyle=1,ystyle=1, $ /overplot,/fill,/isotropic;,POSITION=[ 0.14, 0.08, 0.82, 0.92 ] ;,position=position,zrange=[0.,2.6],zstyle= 1,c_colors = c_colors map_grid,/label,latlab = 175,color=130 map_continents,color=130 ;color bar mxx = MAX(ped2-ped, MIN = mnn ,/nan) * 0.98 delta_lev = ( mxx - mnn ) / FLOAT( nl ) level = FINDGEN( nl ) * delta_lev + mnn ;colorbar array consists of levels previously determined cbar_arr = FLTARR( 2, N_ELEMENTS( level ) ) cbar_arr[0,*] = level cbar_arr[1,*] = level CONTOUR, cbar_arr, [1,2], level, $ NLEVELS = nl, C_COLORS = INDGEN( nl ), XSTYLE = 5, YSTYLE = 5, $ /ZSTYLE, /FILL, BACKGROUND = nl - 1, COLOR = 0, $ POSITION = [ 0.82 + 0.03, 0.03, $ 0.82 + 0.06, 0.945],$ LEVELS = level, /NOERASE, XTICKS = 1, XTICKNAME = [' ', ' ', ' '] AXIS, /YAXIS, charsize=1., charthick=2, /DATA, COLOR = 0, $ /ZSTYLE, /YSTYLE device,/close set_plot,'x' ;openw,33,'gden_400km_1996_M107_3.dat' ;printf,33,'global density at 400km (ng/m3)' ;for n = 0, ntimes-1 do begin ;printf,33,gden(n) ;endfor ;close,33 stop ;stru = {smtime:0.,density:fltarr(72,36,29),tn:fltarr(72,36,29),qjoule:fltarr(72,36,29),$ ; qjoule_integ:fltarr(72,36),sheating:fltarr(72,36,29),zg:fltarr(72,36,29)} ;openw, 111, '2003300_311fism_all.dat' ;printf, 111, 'smtime(day) DEN TN qjoule qjoule_integ sheating zg' ;for i = 0, ntimes-1 do begin ;stru.smtime = smtime(i) ;stru.density = mden(*,*,*,i) ;stru.tn = mtn(*,*,*,i) ;stru.qjoule = mqjoule(*,*,*,i) ;stru.qjoule_integ = mqjoule_integ(*,*,i) ;stru.sheating = msheating(*,*,*,i) ;stru.zg = mzg(*,*,*,i) ;printf, 111, stru ;endfor ;close, 111 int_jh = fltarr(ntimes) & int_jh(*)=0. int_jh2 = fltarr(ntimes) & int_jh2(*)=0. height_jh = fltarr(nlon, nlat,ntimes) & height_jh(*,*,*) =0. ;height integrated joule heating for n=0,ntimes-1 do begin print,n for i = 0, nlon-1 do begin for j=0, nlat-1 do begin done = 0 for k = 1, nlev-2 do begin height_jh(i,j,n) = height_jh(i,j,n)+mqjoule(i,j,k,n)*mden(i,j,k,n)*(mZ(i,j,k+1,n)-mZ(i,j,k-1,n))/2. int_jh(n) = int_jh(n) + $ mqjoule(i,j,k,n)*mden(i,j,k,n)*(mZ(i,j,k+1,n)-mZ(i,j,k-1,n))/2.*cos(lat(j)*!pi/180.) endfor ;k:lev area2 = (Ri+mZ(i,j,nlev-1,n))^2*dlat*dlon int_jh2(n) = int_jh2(n)+mqjoule_integ(i,j,n)*cos(lat(j)*!pi/180.)*area2 endfor ;j:lat endfor ;i:lon endfor ;n:ntimes int_jh2 = int_jh2*1.e-16 ; GW ; Z was in km, and convert from ergs/s*1e5 to GW by a factor of 1.e-11 area = (Ri+mZ(0,0,nlev-1,ntimes-1))^2*dlat*dlon ;YH 06/13/11 int_jh = int_jh*area*1.e-16 ;GW openw, 11, 'hourly_global_JH_1996_M107_300_366.dat' printf, 11, 'smtime(day) int_jh2(GW) int_jh(GW)' for i = 0, ntimes-1 do begin printf, 11, smtime(i), int_jh2(i), int_jh(i) endfor close, 11 ;set_plot,'ps' ;device,filename = 'ijh1996_new.eps' ;plot, smtime_daily,int_jh_daily,ytitle = 'Joule heating(GW)',xtitle='Days',xr=[0,365],xstyle =1 ;device,/close ;set_plot,'X' set_plot,'ps' device, filename ='height_jh_smax.ps' ;device,filename='../tiegcm_plots/JH_smin_2008_F107_70_hp_20_cpcp_50_F_cut_2.ps' device,decomposed=0,/color loadct, 34,ncolors=100,bottom=1 c_colors = indgen(100)+1 position = [0.1,0.1,0.9,0.75] imax = max(height_jh(*,18:35,ntimes-1)) imin = min(height_jh(*,18:35,ntimes-1)) lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180. height_jh2 = fltarr(nlon+1,nlat,ntimes) & height_jh2(0:71,*,*)= height_jh & height_jh2(72,*,*)=height_jh(0,*,*) map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[50.,-180,90.,180] contour, height_jh2(*,18:35,ntimes-1),lon2,lat(18:35),/overplot,/fill,/isotropic,nlevels=100,position=position,zrange=[0.,2.6],zstyle = 1,c_colors = c_colors map_grid,/label,latlab = 175,color=130 map_continents,color=130 ;colorbar,bottom=1,ncolors=50,position=[.2,0.85,0.75,0.9] ;colorbar,bottom=1,divisions=20,ncolors=100,position=[0.2,0.85,0.75,0.9],pscolor=pscolor,title='Joule heating(W/m2)' ;contour, height_jh(*,18:35,ntimes-1),lon,lat(18:35),/overplot,/noerase,nlevels=12 device,/close set_plot,'X' stop ;close,/all end