;+
; NAME:
;	put_spectra
;
; PURPOSE:
;       To write outputs
;
; CATEGORY:
;       utility procedure called by main program.
;
; CALLING SEQUENCE:
;       put_spectra,mode,solar_file,out_spectra,sigma_o,sigma_o2,sigma_n2
;
; INPUTS:
;       mode:        ='dev' for development mode, needs to write 
;                           both tabulated spectra and tabulated 
;                           cross section and branching ratio
;                    ='prod' for production mode, only write 
;                           tabulated spectra
;       solar_file:  (='dir_name/file_name'), original input
;                    spectra file name
;       out_spectra: tabulated spectra 
;
; OUTPUTS:
;       sigma_o:      tabulated cross section and branching ratio 
;                     for atomic oxygen.
;       sigma_o2:     tabulated cross section and branching ratio 
;                     for molecular oxygen.
;       sigma_n2:     tabulated cross section and branching ratio 
;                     for molecular nitrogen.
; COMMON BLOCKS:
;	None.
;
; PROCEDURE:
;
; ROUTINES CALLED:
;       None.
;
; MODIFICATION HISTORY:
;       12/20/02, Liying Qian, Initial Version
;
;+
;

pro put_spectra, mode, iscale, f107, f107a,solar_file,out_spectra,sigma_o,sigma_o2,sigma_n, sigma_n2

n_bins=n_elements(out_spectra(0,*))

if strpos(solar_file,'sc21refw') ne -1 then begin 
   if (iscale eq 0) then begin
     file_name='rflux_'+strcompress(n_bins,/remove_all)+'.dat'
     openw,1,file_name
     printf,1,'Reference spectrum--sc21refw, low solar activity,when read into glow, do flux/10^9'
     printf,1,' wave1(A)', ' wave2',' flux(photon cm^-2 s^-1)', $
  	    ' scale1',' scale2'
     format='$(f7.2,2x,f7.2,2x,e10.3,2(2x,e10.3))'
     for i=0,n_bins-1 do begin
        printf,1,format,out_spectra[0,i],out_spectra[1,i],out_spectra[2,i], $
                 out_spectra[3,i]*1.e-6,out_spectra[4,i]*1.e-6
     endfor
     close,1
   endif else begin
     b1=[1.0, 0.0138, 0.005]
     b2=[1.0, 0.59425, 0.3811]

     r1=b1(0) + b1(1)*(f107a-71.5) + b1(2)*(f107-f107a+3.9)
     r2 =  b2(0) + b2(1)*(f107a-71.5) + b2(2)*(f107-f107a+3.9)

     flux=fltarr(1,n_bins)
     flux[0,*] = out_spectra[2,*] + (r1-1.)*out_spectra[3,*] +(r2-1)*out_spectra[4,*]

     file_name='flux_'+strcompress(n_bins,/remove_all)+'.dat'
     openw,1,file_name
     printf,1,'Scaled spectrum--sc21refw (f107='+strcompress(f107,/remove_all)+', f107a='+strcompress(f107a,/remove_all)+'), when read into glow, do flux/10^9'
     printf,1,' wave1(A)', ' wave2',' flux(photon cm^-2 s^-1)'
     format='$(f7.2,2x,f7.2,2x,e10.3)'
     for i=0,n_bins-1 do begin
        printf,1,format,out_spectra[0,i],out_spectra[1,i],flux[0,i]
     endfor
     close,1
   endelse

endif else if strpos(solar_file,'f79050') ne -1 then begin 
   file_name='xflux_'+strcompress(n_bins,/remove_all)+'.dat'
   openw,1,file_name
   printf,1,'Reference spectrum--f79050, high solar activity,when read into glow, do flux/10^9'
   printf,1,' wave1(A)', ' wave2',' flux(photon cm^-2 s^-1)'
   format='$(f7.2,2x,f7.2,2x,e10.3)'
   for i=0,n_bins-1 do begin
      printf,1,format,out_spectra[0,i],out_spectra[1,i],out_spectra[2,i]
   endfor
   close,1

endif else if strpos(solar_file,'ref_min_27day_11yr') ne -1 then begin 

   file_name='woodsflux_'+strcompress(n_bins,/remove_all)+'.dat'
   openw,1,file_name

   if (iscale eq 0) then begin
     printf,1,'Reference spectrum--ref_min_27day_11yr, when read into glow, do flux/10^9'
     printf,1,' wave1(A)', ' wave2',' flux(photon cm^-2 s^-1)', $
	    ' 27-day variation',' 11-year variation'
     format='$(f7.2,2x,f7.2,2x,e10.3,2(2x,e10.3))'
     for i=0,n_bins-1 do begin
        printf,1,format,out_spectra[0,i],out_spectra[1,i],out_spectra[2,i], $
               out_spectra[3,i],out_spectra[4,i]
     endfor

   endif else begin

     flux=fltarr(1,n_bins)
     flux[0,*]=out_spectra[2,*]+(f107-f107a)*out_spectra[3,*]+f107a*out_spectra[4,*]

     printf,1,'Scaled spectrum--ref_min_27day_11yr(f107='+strcompress(f107,/remove_all)+', f107a='+strcompress(f107a,/remove_all)+'), when read into glow, do flux/10^9'
     printf,1,' wave1(A)', ' wave2',' flux(photon cm^-2 s^-1)'
     format='$(f7.2,2x,f7.2,2x,e10.3)'
     for i=0,n_bins-1 do begin
        printf,1,format,out_spectra[0,i],out_spectra[1,i],flux[0,i]
     endfor
   endelse

   close,1

endif else if ((strpos(solar_file,'.ncdf') ne -1) or (strpos(solar_file,'.nc') ne -1)) then begin 
   read_netcdf,solar_file,data_old
   n_bins=n_elements(out_spectra[0,*])
   n_days=n_elements(out_spectra[*,0])-2
   data_new={DATE:lonarr(n_days),UTTIME:dblarr(n_days),YFRAC:dblarr(n_days),WAVE1:dblarr(n_bins), $
             WAVE2:dblarr(n_bins),SP_FLUX:dblarr(n_bins,n_days)}
   if ( (strpos(solar_file,'L3A') ne -1) or (strpos(solar_file,'FISM') ne -1) ) then begin
;      data_new.UTTIME=double(0.5*(data_old.START_TIME+data_old.STOP_TIME)/(3600.))
      data_new.UTTIME=double(data_old.TIME/3600.)
   endif else begin
      data_new.UTTIME=12.
   endelse
   attr_file='./input/see__L3.attr'
   yfrac=dblarr(n_days)
   for i=0,n_days-1 do begin
      iyear=data_old.DATE[i]/1000
      if ( (iyear-iyear/4*4) eq 0) then begin
         yfrac[i]=double(data_old.DATE[i]/1000)+double(((data_old.DATE[i]-(data_old.DATE[i]/1000)*1000)-1+ $
	                   double(data_new.UTTIME[i]/24.))/366.)
      endif else begin
         yfrac[i]=double(data_old.DATE[i]/1000)+double(((data_old.DATE[i]-(data_old.DATE[i]/1000)*1000)-1+ $
	                   double(data_new.UTTIME[i]/24.))/365.)
      endelse
   endfor
   data_new.YFRAC=yfrac
   data_new.DATE=data_old.DATE
   ; 12/14/2004, I don't understand why I wrote the following commented code
   ; I commented the following code and reverse the data (so that data is
   ; stored iin the order of from longer wavelength to shorter wavelength)
   ; in correct way.
   ;
   ;data_new.WAVE1=reverse(out_spectra[0,*],0)
   ;data_new.WAVE2=reverse(out_spectra[1,*],0)
   ;for i=0,n_days-1 do begin
   ;   data_new.SP_FLUX[*,i]=reverse(out_spectra[i+2,*],0)
   ;endfor
   ;
   temp_array=reverse(out_spectra,2)
   data_new.WAVE1=temp_array[0,*]
   data_new.WAVE2=temp_array[1,*]
   for i=0,n_days-1 do data_new.SP_FLUX[*,i]=temp_array[i+2,*]
   ; end modification on 12/14/2004

   result=strsplit(solar_file,'/',/extract)
   basename=result[n_elements(result)-1]
   result=strsplit(basename,'.',/extract)
   fn=result[0]+'_'+strcompress(n_bins,/remove_all)+'.nc'
   file_name=fn
   print,data_new
   write_netcdf,data_new,file_name,att_file=attr_file

;write Juan's data
endif else if ((strpos(solar_file,'HaoBinsEUV0Fmt.txt') ne -1) ) then begin 
   nrows=2435
   ncolumns=34
   d=dblarr(ncolumns,nrows)
   a=' '
   close,1
   openr,1,solar_file
   readf,1,a
   readf,1,d
   close,1
   n_bins=n_elements(out_spectra[0,*])
   n_days=n_elements(out_spectra[*,0])-2
   data_new={DATE:lonarr(n_days),UTTIME:dblarr(n_days),YFRAC:dblarr(n_days),WAVE1:dblarr(n_bins), $
             WAVE2:dblarr(n_bins),SP_FLUX:dblarr(n_bins,n_days)}
   ;
   data_new.UTTIME=0.
   ;
   year=lonarr(n_days)
   month=lonarr(n_days)
   day=lonarr(n_days)
   year=long(d[1,*]/10000)
   month=long((long(d[1,*])-long(year)*long(10000))/long(100))
   day=(long(d[1,*])-long(year)*long(10000))-month*long(100)
   for m=0,n_days-1 do begin
;      print,year[m],month[m],day[m]
      data_new.date[m]=ymd2yyyydoy(year[m],month[m],day[m])
;      print,data_new.date[m]
   endfor
   ;
   yfrac=dblarr(n_days)
   for i=0,n_days-1 do begin
      iyear=data_new.DATE[i]/1000
      if ( (iyear-iyear/4*4) eq 0) then begin
         yfrac[i]=double(data_new.DATE[i]/1000)+double(((data_new.DATE[i]-(data_new.DATE[i]/1000)*1000)-1+ $
	                   double(data_new.UTTIME[i]/24.))/366.)
      endif else begin
         yfrac[i]=double(data_new.DATE[i]/1000)+double(((data_new.DATE[i]-(data_new.DATE[i]/1000)*1000)-1+ $
	                   double(data_new.UTTIME[i]/24.))/365.)
      endelse
   endfor
   data_new.YFRAC=yfrac
   ; 12/14/2004, I don't understand why I wrote the following commented code
   ; I commented the following code and reverse the data (so that data is
   ; stored iin the order of from longer wavelength to shorter wavelength)
   ; in correct way.
   ;
   ;data_new.WAVE1=reverse(out_spectra[0,*],0)
   ;data_new.WAVE2=reverse(out_spectra[1,*],0)
   ;for i=0,n_days-1 do begin
   ;   data_new.SP_FLUX[*,i]=reverse(out_spectra[i+2,*],0)
   ;endfor
   ;
   temp_array=reverse(out_spectra,2)
   data_new.WAVE1=temp_array[0,*]
   data_new.WAVE2=temp_array[1,*]
   for i=0,n_days-1 do data_new.SP_FLUX[*,i]=temp_array[i+2,*]
   ; end modification on 12/14/2004

   result=strsplit(solar_file,'/',/extract)
   basename=result[n_elements(result)-1]
   result=strsplit(basename,'.',/extract)
   fn=result[0]+'_'+strcompress(n_bins,/remove_all)+'.nc'
   file_name=fn
   attr_file='./input/see__L3.attr'
   write_netcdf,data_new,file_name,att_file=attr_file
endif else begin
   file_name='flux_'+strcompress(n_bins,/remove_all)+'.dat'
   openw,1,file_name
   printf,1,'solar spectra, when read into glow, do flux/10^9'
   printf,1,'wave1(A) ', 'wave2',' fluxes(photon cm^-2 s^-1)'
   result=size(out_spectra)
   n_columns=result[1]
;   fmt='$(2f8.2,12(e11.3,:))'
;   fmt='$(2f8.2,4320(e11.3))'    ; fism 3-day flare spectra
;   fmt='$(2f8.2,28800(e11.3))'    ; fism LWS 2010 20day preflare spectra
;   fmt='$(2f8.2,1500(e11.3))'    ; fism LWS 2010 25 hours control spectra
   fmt='$(2f8.2,5760(e11.3))'    ; fism LWS 2010 4 days after flare
   for i=0,n_bins-1 do begin
      ;printf,1,out_spectra[0,i],out_spectra[1,i],out_spectra[2,i]
      printf,1,fmt,out_spectra[0,i],out_spectra[1,i],out_spectra[2:n_columns-1,i]
   endfor
   close,1
endelse

; Write cross section and branching ratio data

if mode eq 'dev' then begin

  file_name='sigma_'+strcompress(n_bins,/remove_all)+'.dat'
  openw,1,file_name
  printf,1, 'O absorption cross section and ionization branching ratios.'
  printf,1,format = '(a63)','Wavelength Bins(A)     Abs    TotIon      4s      2Do      2Po '
  format="$(f7.2,2x,f7.2,2x,f10.5,4(2x,f7.3))"
  for i=0,n_bins-1 do begin
     printf,1,format, out_spectra[0:1,i],sigma_o(0:4,i)
  endfor

  printf,1, '  '
  printf,1, 'O2 absorption cross section and ionization/dissociation branching ratios.'
  printf,1,format = '(a72)','Wavelength Bins(A)     Abs    TotIon      X     (O2+)*     D.I.    Disso.'
  format="$(f7.2,2x,f7.2,2x,f10.5,5(2x,f7.3))"
  for i=0,n_bins-1 do begin
     printf,1,format, out_spectra[0:1,i],sigma_o2(0:5,i)
  endfor

  printf,1, '  '
  printf,1, 'N2 absorption cross section and ionization/dissociation branching ratios.'
  printf,1,format = '(a73)','Wavelength Bins(A)     Abs    TotIon      X     (N2+)*     D.I.     Disso.'
  format="$(f7.2,2x,f7.2,2x,f10.5,5(2x,f7.3))"
  for i=0,n_bins-1 do begin
     printf,1,format, out_spectra[0:1,i],sigma_n2(0:5,i)
  endfor

  printf,1, '  '
  printf,1, 'N absorption cross section.'
  printf,1,format = '(a26)','Wavelength Bins(A)     Abs'
  format="$(f7.2,2x,f7.2,2x,f10.5)"
  for i=0,n_bins-1 do begin
     printf,1,format, out_spectra[0:1,i],sigma_n(0,i)
  endfor

  close,1

  ; write cross section and branching ratio data in the format as
  ; array constructor for T*GCM)

  temp_o=reverse(sigma_o,2)
  temp_o2=reverse(sigma_o2,2)
  temp_n2=reverse(sigma_n2,2)
  temp_n=reverse(sigma_n,2)

  file_name='sigma_'+strcompress(n_bins,/remove_all)+'_tgcm.dat'
  openw,1,file_name
  n_loop=fix(n_bins/5)
  residual=n_bins-fix(n_bins/5)*5

  printf,1,'      sigeuv(1,:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o2[0,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o2[0,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      sigeuv(2,:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[0,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[0,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      sigeuv(3,:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_n2[0,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_n2[0,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      BPhotonI(1,:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o2[1,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o2[1,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      BPhotonI(2,:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[1,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[1,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      BPhotonI(3,:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_n2[1,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_n2[1,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      brop2pPh(:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[4,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[4,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      brop2dPh(:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[3,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[3,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      brop4sPh(:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[2,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o[2,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      bro2DPh(:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o2[5,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o2[5,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      brn2DPh(:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_n2[5,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_n2[5,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      bro2DIPh(:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_o2[4,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_o2[4,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      brn2DIPh(:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_n2[4,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_n2[4,ind1:ind2]
  endif

  printf,1,' '
  printf,1,'      sigin4s(:) = (/'
  for i=0,n_loop-1 do begin
     ind1=i*5
     ind2=i*5+4
     printf,1,format="('     |         ',5(e9.2,','))",temp_n[0,ind1:ind2]
  endfor
  if (residual ne 0) then begin
     ind1=ind2+1
     ind2=ind2+residual
     printf,1,format="('     |         ',5(e9.2,','))",temp_n[0,ind1:ind2]
  endif

  close,1

  ; write solar flux/cross section for diagnosis 
  ; (used by program deposition.pro to look at ionization rate)

  file_name='comp_'+strcompress(n_bins,/remove_all)+'.dat'
  openw,1,file_name
  printf,1,n_bins
  printf,1,'Output of bin.pro, wave length in A, cross section in megabarns'
  printf,1,'  wavelength ', '  flux  ',' O-I ','  O-A ','  O2-I ',$
  	 '  O2-A ','   N2-I ',' N2-A '
  format='$(f7.2,2x,f7.2,2x,e10.3,6(2x,f9.3))'
  for i=0,n_bins-1 do begin
     printf,1,format,out_spectra[0:2,i],sigma_o[0,i]*sigma_o[1,i], sigma_o[0,i],$
		     sigma_o2[0,i]*sigma_o2[1,i], sigma_o2[0,i],$
		     sigma_n2[0,i]*sigma_n2[1,i],sigma_n2[0,i]
  endfor
  close,1
endif

return
end