source: trunk/mesoscale/PLOT/RESERVE/obsolete/read_dat.pro @ 111

pro read_dat, output, x, y, z, t, desc_field, nodata=nodata


;---------------------------------------------------------------------------
;
; Read the files generated by ARWpost
;
; Please create symbolic link 'input.ctl' and 'input.dat'
; --- or rename files
;
; - output contains all fields
; - x,y,z,t contains coordinates
; - desc_field contains strings describing fields
; - if nodata is not blank, do not read the fields
;
;---------------------------------------------------------------------------
; A. Spiga, August 2007
;---------------------------------------------------------------------------


;
; info for localtime
;
interv=0.
openr, 99, 'timefil'
readf,99,interv
close, 99



; get preliminar info on coordinates
;
wnx=1
wny=1
wnz=1
wnt=1
txt=''
param=[0.,0.,0.,0.,0.,0.,0.,0.]
openr, 1, 'input.ctl'
while not eof(1) do begin
    readf,1,txt
    match=STRPOS(txt,'linear')
        if (match ne -1) then begin
                wnx=(STRPOS(txt,'xdef') eq -1)*wnx  ;0 if linear
                wny=(STRPOS(txt,'ydef') eq -1)*wny  ;0 if linear
                wnz=(STRPOS(txt,'zdef') eq -1)*wnz  ;0 if linear 
                wnt=(STRPOS(txt,'tdef') eq -1)*wnt  ;0 if linear
                paramstart=STRSPLIT(txt, 'r', LENGTH=paramlength)
                paramchar=STRMID(txt,paramstart(1),paramlength(1))
                READS, paramchar, start, step
                param=param+[start*(STRPOS(txt,'xdef') ne -1),$
                        step*(STRPOS(txt,'xdef') ne -1),$
                        start*(STRPOS(txt,'ydef') ne -1),$
                        step*(STRPOS(txt,'ydef') ne -1),$
                        start*(STRPOS(txt,'zdef') ne -1),$
                        step*(STRPOS(txt,'zdef') ne -1),$
                        start*(STRPOS(txt,'tdef') ne -1),$
                        step*(STRPOS(txt,'tdef') ne -1)]
        endif   
endwhile
close, 1
timebegin=param(6)
        
;---------------------------------
; READ .CTL FILE INFORMATIONS
;---------------------------------

info=read_ascii('input.ctl', missing_value=1e37)
infodat=info.field1

;
; second column : dimensions
;

w=where((infodat(0,*) eq 1e37) and (infodat(1,*) lt 1e37))
dimensions=reform(infodat(1,w))

nx=round(dimensions(0)) 
ny=round(dimensions(1)) 
nz=round(dimensions(2)) 
nt=round(dimensions(3)) 
fields=round(dimensions(4)) 
vertical=intarr(fields)
for i=0,fields-1 do begin
        vertical(i)=round(dimensions(5+i))
endfor


;
; first column : coordinates
;

w=where((infodat(1,*) eq 1e37) and (infodat(0,*) lt 1e37))
if (w(0) ne -1) then coordinates=reform(infodat(0,w))

if (wnx ne 0) then begin
        wnx=nx
        x=coordinates(0:nx-1)
endif else begin
        xmin=param(0)
        xmax=param(0)+param(1)*(nx-1) MOD 360
        ;inc=param(0)+param(1)*(nx-1)   ;cas mercator
        ;x=start+inc*findgen(nx) 
        x=xmin+(xmax-xmin)*findgen(nx)/(nx-1)
        ;print, x
        ;stop
endelse
if (wny ne 0) then begin
        wny=ny
        y=coordinates(wnx:wnx+wny-1)
endif else begin
        start=param(2)
        step=param(3)
        y=start+step*findgen(ny)
endelse
if (wnz ne 0) then begin
        wnz=nz
        z=coordinates(wnx+wny:wnx+wny+wnz-1)
endif else begin
        start=param(4)
        step=param(5)
        z=start+step*findgen(nz)
endelse
t=interv*findgen(nt)/3700. + timebegin + mean(x)/15.
        

;
; get info on fields
;
openr, 1, 'input.ctl'
toto=5+wnx+1+wny+1+wnz+2+fields
infofields=strarr(toto)
readf, 1, infofields
close, 1

;
desc_field=transpose(infofields(toto-fields:toto-1))


if (n_elements(nodata) eq 0) then begin
;---------------------------------
; READ .DAT FILE DATA
;---------------------------------
; GrADS views gridded data sets as 5-dimensional arrays 
; varying in longitude, latitude, vertical level, variable, and time.
;---------------------------------

print, '----------------'
print, 'reading data ...'
print, '----------------'
print, 'longitudes: ', min(x), max(x), nx, ' points'
print, 'latitudes: ', min(y), max(y), ny, ' points' 
print, 'altitudes: ', min(z), max(z), nz, ' points'
print, 'local times: ', min(t), max(t), nt, ' points'
print, fields, ' recorded fields :'
print, desc_field
print, '----------------'

        ;;very slow method
        ;;----------------      
        ;data=read_ascii('input.dat', missing_value=1e37)
        ;datadat=data.field1

; see http://www.dfanning.com/tips/ascii_column_data.html       
OPENR, lun, 'input.dat', /GET_LUN

;count=0L
output=fltarr(fields,nx,ny,nz,nt)
dummy=0.

for l=0,nt-1 do begin                                           ;time loop
print, 'time loop ...', l+1, ' / ', string(nt,'(I0)')   
        for f=0, fields-1 do begin ;fields (whether 2D or 3D)   ;variable loop
        nzf=vertical(f)
        print, 'reading field ...', f+1, ' / ', string(fields,'(I0)') 
                for k=0,nzf-1 do begin                          ;vertical loop
                        for j=0,ny-1 do begin                   ;latitude loop
                        for i=0,nx-1 do begin                   ;longitude loop

                                ;;very slow method
                                ;output(f,i,j,k,l)=datadat(count)
                                ;count=count+1

                                READF, lun, dummy
                                output(f,i,j,k,l)=dummy

                        endfor
                        endfor
                endfor
        endfor
endfor

; missing values
w=where(output eq 1e20)
if (w(0) ne -1) then begin
        print, 'change missing values ...'
        output[w]=1e37
endif


CLOSE, lun
;SPAWN, '\rm input.dat_tmp'
print, '...done !'
help, /memory
print, '---------'

;if (count ne n_elements(datadat)) then begin
;       print, 'some data were not read !'
;       stop
;endif

;;
;; Manage missing data so as IDL recognize them
;;
;w=where(output eq 1e37)
;if (w(0) ne -1) then output[w]=!Values.F_NAN

endif

end