;---------------------------------------------------
;                 XIAOHUA FANG                      |
;   Laboratory for Atmospheric and Space Physics    |
;         University of Colorado, Boulder           |
;         xiaohua.fang@lasp.colorado.edu            |
;                  11/10/2011                       |
;                STEPHEN BOUGHER
;                 AOSS/ Michigan
;                  11/16/11
;---------------------------------------------------

pro arrays_swb_ncdf3
; --------------------------------------------------------------------------------------------------
    delta =5.         ; grid resolution [degree]: the same in the longitude and latitude directions
    Ntheta=36         ; theta: co-latitude, measured from the North Pole 
    nlat = 36
    Nphi  =72         ; phi: "longitude", measured from the nightside (LT=0)
    nlon = 72
;   filename='impact_case1.dat'
;   filename='impact_case1_new.dat'
    filename='impact_case7.dat'
; --------------------------------------------------------------------------------------------------

    ; constants
    ImpactRate_ = 0
    AveEnergy_  = 1
    
    ;\
    ; read data file
    ;/

    ; array/variable definitions
    str='temporary string'
    buff=fltarr(4)    ; buffer 
    data=fltarr(3, Ntheta, Nphi)
    ; MTGCM array/variable definitions
    ; -- fast variable: longitude  (0 degress at midnight, by 5 degree interavls)
    ; -- slow variable: latitude   (-87.5 to +87.5, by 5 degree intervals)
    opflux=fltarr(72,36)
    opave=fltarr(72,36)
    glat = fltarr(36)
    glon = fltarr(72)
    
    ; read data file
    openr, 1, filename
    for i=1, 6 do begin
        ; skip the heading text
        readf, 1, str
    endfor
    for jP=0, Nphi-1 do begin
        longitude=jP*delta                              ; longitude, measured from the nightside, LT=0 [degree]
        glon(jP) = longitude
        for iT=0, Ntheta-1 do begin
            theta=(iT+0.5)*delta                        ; theta: co-latitude, measured from the North Pole
            latitude=90.-theta                          ; latitude [degree]
            glat(iT)= -90. + theta

            readf, 1, buff

            ; self test only
            if (abs(buff(0)-longitude) gt 0.01) or (abs(buff(1)-latitude) gt 0.01) then begin
                print, '### ERROR ###'
                STOP
            endif  
          
            data(ImpactRate_, iT, jP)=buff(2)           ; O+ impact rate at the exobase altitude [cm^-2 s^-1]
            data(AveEnergy_ , iT, jP)=buff(3)           ; average energy of impacting O+ [eV]
;           opflux(jP,35-iT) = data(ImpactRate_, iT, jP) ; O+ impact rate into MTGCM grid array
;           opave(jP,35-iT) = data(AveEnergy_ , iT, jP)  ; average energy into MTGCM grid array

            if (jP le 35) then begin
                opflux(jP+36,35-iT) = data(ImpactRate_, iT, jP) ; rotate 12-hours; find why!
                opave(jP+36,35-iT) = data(AveEnergy_ , iT, jP)  ; rotate 12-hours; find why!
            endif
            if (jP gt 35) then begin
                opflux(jP-36,35-iT) = data(ImpactRate_, iT, jP) ; rotate 12-hours; find why!
                opave(jP-36,35-iT) = data(AveEnergy_ , iT, jP)  ; rotate 12-hours; find why!
            endif

;            print, longitude, latitude, data(ImpactRate_, iT, jP), data(AveEnergy_ , iT, jP),   $
;                   format='(F5.1, F7.1, 2E12.3)'
        endfor
    endfor
    close, 1
; -------------------------------------------------------------------------------------
;  Print to files new MTGCM arrays
;   output2 = 'OPFLUX.asc'
;   output3 = 'OPAVE.asc'
;  openw, 2, output2
;  for iT=0, Ntheta-1 do for jP=0,Nphi-1 do printf,2,opflux(jP,iT), format='(E12.3)'
;   for iT=0, Ntheta-1 do begin
;       for jP=0, Nphi-1 do begin
;           print,opflux(jp,iT), format='(6E12.3)' 
;       endfor
;   endfor
;  openw, 3, output3
;  for iT=0, Ntheta-1 do for jP=0,Nphi-1 do printf,3,opave(jP,iT), format='(E12.3)'
;   for iT=0, Ntheta-1 do begin
;       for jP=0, Nphi-1 do begin
;           print,opave(jp,iT) 
;       endfor
;   endfor
;   close, 2
;   close, 3
; -------------------------------------------------------------------------------------
; Foster: I'm attaching a procedure that writes a netcdf file with IMF data,
; but the basic idea is, you create a new dataset (nc file), define
; dimensions, variables and attributes, go out of define mode and into
; data mode, and give values to dims and vars.
; ---------------------------------------------------------------------
  ncid = ncdf_create('file.nc',/clobber) ; creat file (overwrite preexist)
 
  latdim = ncdf_dimdef(ncid,'lat',nlat) ; create latitude dimension
  londim = ncdf_dimdef(ncid,'lon',nlon) ; create longitude dimension
 
; Define coord vars (they have same name as corresponding dimensions)
  idlat = ncdf_vardef(ncid,'lat',[latdim],/FLOAT)
  idlon = ncdf_vardef(ncid,'lon',[londim],/FLOAT)
;
; Create a 2d var (lon,lat), with couple of attributes:
; (create as many vars as you want)

; OPFLUX Variable
  id_var1 = ncdf_vardef(ncid,'var1',[londim,latdim],/FLOAT) ;create 2dvar
  ncdf_attput,ncid,id_var1,'long_name','OP Flux7' ; give it long name
  ncdf_attput,ncid,id_var1,'units','cm-2sec-1' ; units attribute

; OPaverage Energy Variable
  id_var2 = ncdf_vardef(ncid,'var2',[londim,latdim],/FLOAT) ;create 2dvar
  ncdf_attput,ncid,id_var2,'long_name','OP Average Energy7' ; give it long name
  ncdf_attput,ncid,id_var2,'units','eV' ; units attribute
 
  ncdf_control,ncid,/endef ; take out of define mode, into data mode
 
  ncdf_varput,ncid,idlat,glat ; give values to latitude coord var
  ncdf_varput,ncid,idlon,glon ; give values to longitude coord var

  ncdf_varput,ncid,id_var1,opflux ; give values to 2d var
  result1 = ncdf_varinq(ncid,id_var1)
  help, result1, /Structure
  ncdf_varput,ncid,id_var2,opave  ; give values to 2d var
  result2 = ncdf_varinq(ncid,id_var2)
  help, result2, /Structure

  ncdf_close,ncid ; close file

; ---------------------------------------------------------------------

 end