function fslt, slt,ut,glon,iflag ; ; If iflag = 1 -> return slt, given ut and glon ; If iflag = 2 -> return ut, given slt and glon ; If iflag = 3 -> return glon, given slt and ut ; (slt and ut always in decimal hours, glon in -180 -> 180) ; case iflag of 1: begin fslt = ut + glon/15. if (fslt lt 0.) then fslt = fslt+24. if (fslt ge 24.) then fslt = fslt-24. end 2: begin fslt = slt - glon/15. if (fslt lt 0.) then fslt = fslt+24. if (fslt gt 24.) then fslt = fslt-24. end 3: begin fslt = (slt-ut) * 15. if (fslt ge 180.) then fslt = fslt-360. if (fslt lt -180.) then fslt = fslt+360. end else: begin print,'>>> fslt: bad iflag=',iflag fslt = 0. end endcase return,fslt end PRO meped_ring_fit ; fit NOAA MEPED electron and proton mean energy and energy flux dayofmon=[31,28,31,30,31,30,31,31,30,31,30,31] linein = '' dummy = fltarr(15) caption = strarr(2,7,7) missing = 0. min_ekev = 30. max_ekev = 1000. nlon = 180 & nlat = 60 ; reassign to TGCM grids nnx=289 & nny=48 min_lat = 30. flux = fltarr(2,7,nlon,nlat) &flux0 = flux glat = fltarr(nlat) for i=0, nlat-1 do glat(i) = 30.5 + i*1. glon = fltarr(nlon) for i=0, nlon-1 do glon(i) = 1. + i*2. gglon = glon for i=0, nlon-1 do if (glon(i) gt 180.) then gglon(i)=glon(i)-360. ext = 'oct*_beltdifffluxes' ;read, 'Enter files to process (or oct*fluxes): ',ext list = findfile(ext,count=nfiles) if n_elements(list) gt 0 then for ll=0,n_elements(list)-1 do begin filein = list[ll] close,1 print,'Reading file: ',filein openr,1,filein ; flux is the number flux for electron and proton at 7 differential energy ; channels over 1 (glat) by 2 (glon) above |30|-degree ; the first index refers to the hemisphere, 1=northern, 2=southern ; the second index refers to the energy band ; 1= 800 to 2500 keV protons (ions) ; 2= 240 to 800 keV protons (ions) ; 3= 80 to 240 keV protons (ions) ; 4= 30 to 80 keV protons (ions) ; 5= 300 to 1000 keV electrons ; 6= 100 to 300 keV electrons ; 7= 30 to 100 keV electrons ; the third index refers to a geographic latitude ranging ; from 30.0 to 90.0. If k is the third index, then the ; center of the latitude box is given by ; latitude=30.5+(k-1)*1. ; it is understood that in the southern hemisphere, the latitude ; is negative ; the fourth index refers to a geographic longitude ranging from ; 0 through 360. If l is the fourth index, then the center of ; the longitude box is given by ; longitude=1.0+(l-1)*2. ; this loop handles the 7 particle species and energy bands for j=0,6 do begin ; this loop handles first the north and then the south for i=0,1 do begin ; read past the 7 captions for m=0,6 do begin readf,1,linein caption(i,j,m) = linein endfor ;print,'caption = ',caption(i,j,4) ; now read the data. The 60 latitude bins are stored in 4 lines of ; 15 entries each - starting with 30-31 and ending with 89-90. The ; first 4 lines are appropriate to longitudes 0 to 2. The 4 line ; format repeats 180 times for the 180 longitude bins ; this loop handles the 180 longitude bins for k=0,nlon-1 do begin for n=0,3 do begin n1=n*15+0 n2=n1+14 readf,1,dummy flux0(i,j,k,n1:n2) = dummy(0:14) endfor ; now convert the data back to actual flux units in particles per cm2 ; per second per steradian. for n=0,59 do begin flux(i,j,k,n)=10.^flux0(i,j,k,n) if (flux(i,j,k,n) lt 10.) then flux(i,j,k,n) = missing endfor ; if a flux is 0.5, it means that the results of the fitting routine ; returned a value of less than or equal to 0.5. ; if a flux is 0.05, it means that there was no data in that latitude- ; longitude interval and the flux was set to a very low default value. ; this ends the do k=1,180 loop endfor ; this ends the do i=1,2 loop endfor ; this ends the do j=1,7 loop endfor N_flux_e = fltarr(2,180,60) & N_flux_p = N_flux_e E_flux_e = fltarr(2,180,60) & E_flux_p = E_flux_e E_mean_e = fltarr(2,180,60) & E_mean_p = E_mean_e for lat=0,nlat-1 do for lon=0,nlon-1 do begin ; for northern hemisphere electrons N_int_30kev = flux(0,6,lon,lat)+flux(0,5,lon,lat)+flux(0,4,lon,lat) N_int_100kev = flux(0,5,lon,lat)+flux(0,4,lon,lat) if (N_int_30kev lt 100. or N_int_100kev eq N_int_30kev) then begin E_flux_e(0,lon,lat) = missing N_flux_e(0,lon,lat) = missing E_mean_e(0,lon,lat) = missing E_mean_e(0,lon,lat) = min_ekev endif else begin ; Finding out the parameters for exponantional fitting ; with f(E) = b(1)*exp(-E*b(2)) b = fltarr(3) b(2) = alog(N_int_30kev/N_int_100kev)/70. if (b(2) eq 0.) then stop if (1./b(2) lt min_ekev) then b(2)=1./min_ekev if (1./b(2) gt max_ekev) then b(2)=1./max_ekev b(1) = b(2)*flux(0,6,lon,lat)/(exp(-30*b(2))-exp(-100.*b(2))) N_flux_e(0,lon,lat) = N_int_30kev E_mean_e(0,lon,lat) = 1./b(2) E_flux_e(0,lon,lat) = N_flux_e(0,lon,lat)*E_mean_e(0,lon,lat)*1.602e-9 ;convert from keV/cm^2/s to mW/m^2/s endelse ; for northern hemisphere proton N_int_30kev = flux(0,3,lon,lat)+flux(0,2,lon,lat)+ $ flux(0,1,lon,lat)+flux(0,0,lon,lat) N_int_80kev = flux(0,2,lon,lat)+flux(0,1,lon,lat)+flux(0,0,lon,lat) N_int_240kev = flux(0,1,lon,lat)+flux(0,0,lon,lat) if (N_int_30kev lt 100. or N_int_80kev eq N_int_30kev) then begin E_flux_p(0,lon,lat) = missing N_flux_p(0,lon,lat) = missing E_mean_p(0,lon,lat) = missing E_mean_p(0,lon,lat) = min_ekev endif else begin ; Finding out the parameters for exponantional fitting ; with f(E) = b(1)*exp(-E*b(2)) b = fltarr(3) b(2) = alog(N_int_30kev/N_int_80kev)/50. if (b(2) eq 0.) then stop if (1./b(2) lt min_ekev) then b(2)=1./min_ekev if (1./b(2) gt max_ekev) then b(2)=1./max_ekev b(1) = b(2)*flux(0,3,lon,lat)/(exp(-30*b(2))-exp(-240.*b(2))) N_flux_p(0,lon,lat) = N_int_30kev E_mean_p(0,lon,lat) = 1./b(2) ;E_flux_p(0,lon,lat) = b(1)/b(2)^2*1.6e-9 E_flux_p(0,lon,lat) = N_flux_p(0,lon,lat)*E_mean_p(0,lon,lat)*1.602e-9 endelse ; for southern hemisphere N_int_30kev = flux(1,6,lon,lat)+flux(1,5,lon,lat)+flux(1,4,lon,lat) N_int_100kev = flux(1,5,lon,lat)+flux(1,4,lon,lat) N_flux_e(1,lon,lat) = N_int_30kev if (N_int_30kev lt 100. or N_int_100kev eq N_int_30kev) then begin E_flux_e(1,lon,lat) = missing N_flux_e(1,lon,lat) = missing E_mean_e(1,lon,lat) = missing E_mean_e(1,lon,lat) = min_ekev endif else begin ; Finding out the parameters for exponantional fitting ; with f(E) = b(1)*exp(-E*b(2)) b = fltarr(3) b(2) = alog(N_int_30kev/N_int_100kev)/70. if (b(2) eq 0) then stop if (1./b(2) lt min_ekev) then b(2) = 1./min_ekev if (1./b(2) gt max_ekev) then b(2)=1./max_ekev b(1) = b(2)*flux(1,6,lon,lat)/(exp(-30*b(2))-exp(-100.*b(2))) N_flux_e(1,lon,lat) = N_int_30kev E_mean_e(1,lon,lat) = 1./b(2) ;E_flux_e(1,lon,lat) = b(1)/b(2)^2*1.6e-9 E_flux_e(1,lon,lat) = N_flux_e(1,lon,lat)*E_mean_e(1,lon,lat)*1.602e-9 endelse ; for southern hemisphere proton N_int_30kev = flux(1,3,lon,lat)+flux(1,2,lon,lat)+ $ flux(1,1,lon,lat)+flux(1,0,lon,lat) N_int_80kev = flux(1,2,lon,lat)+flux(1,1,lon,lat)+flux(1,0,lon,lat) N_int_240kev = flux(1,1,lon,lat)+flux(1,0,lon,lat) if (N_int_30kev lt 100. or N_int_80kev eq N_int_30kev) then begin E_flux_p(1,lon,lat) = missing N_flux_p(1,lon,lat) = missing E_mean_p(1,lon,lat) = missing E_mean_p(1,lon,lat) = min_ekev endif else begin ; Finding out the parameters for exponantional fitting ; with f(E) = b(1)*exp(-E*b(2)) b = fltarr(3) b(2) = alog(N_int_30kev/N_int_80kev)/50. if (b(2) eq 0) then stop if (1./b(2) lt min_ekev) then b(2) = 1./min_ekev if (1./b(2) gt max_ekev) then b(2)=1./max_ekev b(1) = b(2)*flux(1,3,lon,lat)/(exp(-30*b(2))-exp(-240.*b(2))) N_flux_p(1,lon,lat) = N_int_30kev E_mean_p(1,lon,lat) = 1./b(2) E_flux_p(1,lon,lat) = N_flux_p(1,lon,lat)*E_mean_p(1,lon,lat)*1.602e-9 endelse endfor if (ll eq 0) then begin ntime = n_elements(list) ut = fltarr(ntime) & year=ut & jday=ut & month=ut & day=ut e_flux_nh = fltarr(ntime,nlon,nlat) e_ekev_nh = fltarr(ntime,nlon,nlat) e_flux_sh = fltarr(ntime,nlon,nlat) e_ekev_sh = fltarr(ntime,nlon,nlat) p_flux_nh = fltarr(ntime,nlon,nlat) p_ekev_nh = fltarr(ntime,nlon,nlat) p_flux_sh = fltarr(ntime,nlon,nlat) p_ekev_sh = fltarr(ntime,nlon,nlat) data_tgcm = fltarr(ntime,8,nnx,nny) endif e_flux_nh(ll,*,*) = E_flux_e(0,*,*) e_ekev_nh(ll,*,*) = E_mean_e(0,*,*) res = where(e_flux_nh(ll,*,*) ne missing, knt) eflux=reform(e_flux_nh(ll,*,*)) emean=reform(e_ekev_nh(ll,*,*)) if (knt gt 1) then print,'Good data point in e_flux_nh,mean_flux,median_E = ',knt,mean(eflux(res)),median(emean(res)) e_flux_sh(ll,*,*) = E_flux_e(1,*,*) e_ekev_sh(ll,*,*) = E_mean_e(1,*,*) res = where(e_flux_sh(ll,*,*) ne missing, knt) eflux=reform(e_flux_sh(ll,*,*)) emean=reform(e_ekev_sh(ll,*,*)) if (knt gt 1) then print,'Good data point in e_flux_sh,mean_flux,median_E = ',knt,mean(eflux(res)),median(emean(res)) p_flux_nh(ll,*,*) = E_flux_p(0,*,*) p_ekev_nh(ll,*,*) = E_mean_p(0,*,*) res = where(p_flux_nh(ll,*,*) ne missing, knt) pflux=reform(p_flux_nh(ll,*,*)) emean=reform(p_ekev_nh(ll,*,*)) if (knt gt 1) then print,'Good data point in p_flux_nh,mean_flux,median_E = ',knt,mean(eflux(res)),median(emean(res)) p_flux_sh(ll,*,*) = E_flux_p(1,*,*) p_ekev_sh(ll,*,*) = E_mean_p(1,*,*) res = where(p_flux_sh(ll,*,*) ne missing, knt) pflux=reform(p_flux_sh(ll,*,*)) emean=reform(p_ekev_sh(ll,*,*)) if (knt gt 1) then print,'Good data point in p_flux_sh,mean_flux,median_E = ',knt,mean(eflux(res)),median(emean(res)) endfor ;Saving netCDF output id_save = 1 if (id_save ge 1) then begin print,'Saving netCDF output >>>' ntime = ll ncid = ncdf_create('med_ring_oct03.nc',/clobber) time_dim = ncdf_dimdef(ncid,'ntime',ntime) lat_dim = ncdf_dimdef(ncid,'nlat',nlat) lon_dim = ncdf_dimdef(ncid,'nlon',nlon) time_var = ncdf_vardef(ncid,'ut',time_dim,/float) lat_var = ncdf_vardef(ncid,'glat',lat_dim,/float) lon_var = ncdf_vardef(ncid,'glon',lon_dim,/float) ;1-D fields year_var = ncdf_vardef(ncid,'year',time_dim,/short) jday_var = ncdf_vardef(ncid,'jday',time_dim,/short) month_var = ncdf_vardef(ncid,'month',time_dim,/short) day_var = ncdf_vardef(ncid,'day',time_dim,/short) ;3-D fields e_flux_nh_var = ncdf_vardef(ncid,'e_flux_nh',[time_dim,lon_dim,lat_dim],/float) e_ekev_nh_var = ncdf_vardef(ncid,'e_ekev_nh',[time_dim,lon_dim,lat_dim],/float) e_flux_sh_var = ncdf_vardef(ncid,'e_flux_sh',[time_dim,lon_dim,lat_dim],/float) e_ekev_sh_var = ncdf_vardef(ncid,'e_ekev_sh',[time_dim,lon_dim,lat_dim],/float) p_flux_nh_var = ncdf_vardef(ncid,'p_flux_nh',[time_dim,lon_dim,lat_dim],/float) p_ekev_nh_var = ncdf_vardef(ncid,'p_ekev_nh',[time_dim,lon_dim,lat_dim],/float) p_flux_sh_var = ncdf_vardef(ncid,'p_flux_sh',[time_dim,lon_dim,lat_dim],/float) p_ekev_sh_var = ncdf_vardef(ncid,'p_ekev_sh',[time_dim,lon_dim,lat_dim],/float) NCDF_ATTPUT,ncid,lat_var,'standard_name','geographic latitude' NCDF_ATTPUT,ncid,lat_var,'units','degrees' NCDF_ATTPUT,ncid,lon_var,'standard_name','geographico longitude' NCDF_ATTPUT,ncid,lon_var,'units','degrees' NCDF_ATTPUT,ncid,time_var,'standard_name','ut' NCDF_ATTPUT,ncid,time_var,'units','hour' NCDF_ATTPUT,ncid,year_var,'standard_name','year' NCDF_ATTPUT,ncid,year_var,'units',' ' NCDF_ATTPUT,ncid,jday_var,'standard name','Julian Day' NCDF_ATTPUT,ncid,jday_var,'units',' ' NCDF_ATTPUT,ncid,month_var,'standard_name','month' NCDF_ATTPUT,ncid,month_var,'units',' ' NCDF_ATTPUT,ncid,day_var,'standard_name','day' NCDF_ATTPUT,ncid,day_var,'units',' ' NCDF_ATTPUT,ncid,e_flux_nh_var,'long_name','MEPED NH electron energy flux' NCDF_ATTPUT,ncid,e_flux_nh_var,'units','mW/m^2' NCDF_ATTPUT,ncid,e_ekev_nh_var,'long_name','MEPED NH electron mean energy' NCDF_ATTPUT,ncid,e_ekev_nh_var,'units','keV' NCDF_ATTPUT,ncid,e_flux_sh_var,'long_name','MEPED SH electron energy flux' NCDF_ATTPUT,ncid,e_flux_sh_var,'units','mW/m^2' NCDF_ATTPUT,ncid,e_ekev_sh_var,'long_name','MEPED SH electron mean energy' NCDF_ATTPUT,ncid,e_ekev_sh_var,'units','keV' NCDF_ATTPUT,ncid,p_flux_nh_var,'long_name','MEPED NH proton energy flux' NCDF_ATTPUT,ncid,p_flux_nh_var,'units','mW/m^2' NCDF_ATTPUT,ncid,p_ekev_nh_var,'long_name','MEPED NH proton mean energy' NCDF_ATTPUT,ncid,p_ekev_nh_var,'units','keV' NCDF_ATTPUT,ncid,p_flux_sh_var,'long_name','MEPED SH proton energy flux' NCDF_ATTPUT,ncid,p_flux_sh_var,'units','mW/m^2' NCDF_ATTPUT,ncid,p_ekev_sh_var,'long_name','MEPED SH proton mean energy' NCDF_ATTPUT,ncid,p_ekev_sh_var,'units','keV' ; Add global attributes: NCDF_ATTPUT, ncid, /GLOBAL,'title','MEPED Data File = '+ext NCDF_CONTROL,ncid,/ENDEF ; Write data output NCDF_VARPUT,ncid,lat_var,glat ;NCDF_VARPUT,ncid,lon_var,gglon NCDF_VARPUT,ncid,lon_var,glon NCDF_VARPUT,ncid,time_var,ut NCDF_VARPUT,ncid,year_var,year NCDF_VARPUT,ncid,jday_var,jday NCDF_VARPUT,ncid,month_var,month NCDF_VARPUT,ncid,day_var,day NCDF_VARPUT,ncid,e_flux_nh_var,e_flux_nh NCDF_VARPUT,ncid,e_ekev_nh_var,e_ekev_nh NCDF_VARPUT,ncid,e_flux_sh_var,e_flux_sh NCDF_VARPUT,ncid,e_ekev_sh_var,e_ekev_sh NCDF_VARPUT,ncid,p_flux_nh_var,p_flux_nh NCDF_VARPUT,ncid,p_ekev_nh_var,p_ekev_nh NCDF_VARPUT,ncid,p_flux_sh_var,p_flux_sh NCDF_VARPUT,ncid,p_ekev_sh_var,p_ekev_sh NCDF_CLOSE,ncid endif STOP end