;---------------------------------------------------
;                STEPHEN BOUGHER
;                 AOSS/ Michigan
;                  11/08/12
; 4-arrays to readin (Z, T, U, V): 2-D --> X(zimx,zjmx) = X(72,36)
; Combine into one larger array: 3-D --->  X(72,36,4)
;---------------------------------------------------

pro arrays_swb_vtgcmlbc_ncdf
; --------------------------------------------------------------------------------------------------
    delta =5.         ; grid resolution [degree]: the same in the longitude and latitude directions
    Ntheta=36         ; theta: latitude, measured from the South Pole 
    nlat = 36
    Nphi  =72         ; phi: "longitude", measured from the nightside (LT=0)
    nlon = 72
    filename1='GeoHtout.txt'
    filename2='Tempout.txt'
;   filename3='Zonalout.txt'
    filename3='Zonalout.new.txt'
    filename4='Meriout.txt'
; --------------------------------------------------------------------------------------------------

    ;\
    ; read data files
    ;/

    ; array/variable definitions
    str='temporary string'
;   buff=fltarr(1)             ; buffer 
    data=fltarr(Nphi,Ntheta,4) ; destination array to be assembled
    ; VTGCM 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)
    glat = fltarr(36)
    glon = fltarr(72)
    
    ; read data files
    openr, 1, filename1
    openr, 2, filename2
    openr, 3, filename3
    openr, 4, filename4
;   for i=1, 6 do begin
;       ; skip the heading text
;       readf, 1, str
;   endfor
    for iT=0, Ntheta-1 do begin
            theta=(iT+0.5)*delta                        ; theta: latitude, measured from the South Pole
            latitude= -90.+theta                        ; latitude [degree]
            glat(iT)=  latitude
        for jP=0, Nphi-1 do begin
            longitude=jP*delta                          ; longitude, measured from the nightside, LT=0 [degree]
            glon(jP) = longitude

            readf, 1, buff
            data(jP, iT, 0)=buff                       ; Z-geohts (meters)
            readf, 2, buff
            data(jP, iT, 1)=buff                       ; TN-temperatures (K)
            readf, 3, buff
            data(jP, iT, 2)=buff                       ; UN- zonal winds  (m/sec)
            readf, 4, buff
            data(jP, iT, 3)=buff                       ; VN- meridional winds  (m/sec)

             print, longitude, latitude, data(jP ,iT, 0), data(jP ,iT, 1),data(jP ,iT, 2),$
	            data(jP ,iT, 3),format='(6F10.1)'
        endfor
    endfor
    close, 1
    close, 2
    close, 3
    close, 4
; -------------------------------------------------------------------------------------
;  Print to files new VTGCM assembled array data(72, 36, 4)
    output1 = 'data4.asc'
    openw, 5, output1
    for iT=0, Ntheta-1 do begin
        for jP=0, Nphi-1 do begin
             print, data(jP ,iT, 0), data(jP ,iT, 1),data(jP ,iT, 2),data(jP ,iT, 3),$
                    format='(4F10.1)'
        endfor
    endfor
    close, 5
; -------------------------------------------------------------------------------------
; 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