;
;  Read in the ecmwf data
;
pro get_ecmwf_data,site,searchdate
;
;  common containing the data
;
common ecmwf,newpress,newuwind,newvwind,newtemp,newshum,newliq,newice,newcfrac,$
  newrhum,newomega,numtimes,base_time,time_offset,numheights,newheight
;
;  Date to be analyzed
;
curdate=strcompress(string(searchdate),/remove_all)
month=strmid(curdate,4,2)
site=strcompress(site,/remove_all)
;
;  Open the netcdf file
;
if site eq 'sgp' then begin
  armdatadir='/data/mace3/arm_data/sgp/'
  datastream='sgpecmwf/'
  arm_file='sgpecmwfvarX1.c1.'+curdate+'*cd*'
endif else if site eq 'nsa' then begin
  armdatadir='/data/mace2/arm_data/nsa/'
  datastream='nsaecmwf/'
  arm_file='nsaecmwfvarX1.c1.'+curdate+'*cd*'
endif else if site eq 'twpc1' then begin
  armdatadir='/data/mace/arm_data/twp/Manus_C1/'
  datastream='twpecmwfman/'
  arm_file='twpecmwfmanvarX1.c1.'+curdate+'*cd*'
endif else if site eq 'twpc2' then begin
  armdatadir='/data/mace/arm_data/twp/Nauru_C2/'
  datastream='twpecmwfnau/'
  arm_file='twpecmwfnauvarX1.c1.'+curdate+'*cd*'
endif
;
;  Create and change to the data directory
;
path=armdatadir+datastream
pushd,path
;
;  Unzip the data files
;
unix_command='gunzip -fq '+arm_file
spawn,unix_command
;
;  Get a list of the matching file names
;
print, 'filename ',arm_file
files=findfile(arm_file,count=num_files)
;
;  Loop through all the matching files in the directory
;
for i=0,num_files-1 do begin
;
;  Don't try to process any gif files that may be in the directory
;
  if strpos(files[i],'gif') eq -1 and $
    strpos(files[i],'png') eq -1 then begin 
    print,i,files[i]
;
;  Open the netcdf file
;
    cdfid=ncdf_open(files[i])
;
;  Get the dimension id's
;
    time_did=ncdf_dimid(cdfid,'time')
    levels_did=ncdf_dimid(cdfid,'levels')
    domains_did=ncdf_dimid(cdfid,'domains')
;
;  Get the dimensions
;
    if time_did ge 0 then ncdf_diminq,cdfid,time_did,charstring,numtimes
    if levels_did ge 0 then ncdf_diminq,cdfid,levels_did,charstring,numlevels
    if domains_did ge 0 then ncdf_diminq,cdfid,domains_did,charstring,numdomains
;
;  Get the variable id's
;
    base_time_id=ncdf_varid(cdfid,'base_time')  
    time_offset_id=ncdf_varid(cdfid,'time_offset')  ;time(time)
    domains_id=ncdf_varid(cdfid,'domains')  ;domains(domains)
    lev_id=ncdf_varid(cdfid,'lev')  ;lev(levels)
    idat_id=ncdf_varid(cdfid,'idat')  ;initial date of the forcast(time)
    itim_id=ncdf_varid(cdfid,'itim')  ;initial time of the forcast(time)
    vdat_id=ncdf_varid(cdfid,'vdat')  ;verifying date of the forcast(time)
    vtim_id=ncdf_varid(cdfid,'vtim')  ;verifying time of the forcast(time)
    press_id=ncdf_varid(cdfid,'p')  ;pressure(time,levels,domains)
    uwind_id=ncdf_varid(cdfid,'u')  ;zonal wind component(time,levels,domains)
    ;meridional wind component(time,levels,domains)
    vwind_id=ncdf_varid(cdfid,'v')  
    temp_id=ncdf_varid(cdfid,'T')  ;temperature(time,levels,domains)
    shum_id=ncdf_varid(cdfid,'q')  ;specific humidity(time,levels,domains)
    liq_id=ncdf_varid(cdfid,'l')  ;specific cloud liquid water vapor content(")
    ice_id=ncdf_varid(cdfid,'i')  ;specific cloud ice content(")
    cfrac_id=ncdf_varid(cdfid,'a')  ;cloud fraction(")
    rhum_id=ncdf_varid(cdfid,'R')  ;relative humidiy(")
    omega_id=ncdf_varid(cdfid,'w')  ;vertical velocity(")
;
;  Get the variables
;
    ;base time
    if base_time_id ge 0 then ncdf_varget,cdfid,base_time_id,base_time  
    ;time offset
    if time_offset_id ge 0 then ncdf_varget,cdfid,time_offset_id,time_offset  
    if domains_id ge 0 then ncdf_varget,cdfid,domains_id,domains  
    if lev_id ge 0 then ncdf_varget,cdfid,lev_id,lev  ;flux levels
    if idat_id ge 0 then ncdf_varget,cdfid,idat_id,idat  ;initial date
    if itim_id ge 0 then ncdf_varget,cdfid,itim_id,itim  ;forecast time
    if vdat_id ge 0 then ncdf_varget,cdfid,vdat_id,vdat  ;verify date
    if vtim_id ge 0 then ncdf_varget,cdfid,vtim_id,vtim  ;verify time
    ;domain,level,time
    if press_id ge 0 then ncdf_varget,cdfid,press_id,press 	
    if uwind_id ge 0 then ncdf_varget,cdfid,uwind_id,uwind  
    if vwind_id ge 0 then ncdf_varget,cdfid,vwind_id,vwind
    if temp_id ge 0 then ncdf_varget,cdfid,temp_id,temp
    if shum_id ge 0 then ncdf_varget,cdfid,shum_id,shum
    if liq_id ge 0 then ncdf_varget,cdfid,liq_id,liq
    if ice_id ge 0 then ncdf_varget,cdfid,ice_id,ice
    if cfrac_id ge 0 then ncdf_varget,cdfid,cfrac_id,cfrac
    if rhum_id ge 0 then ncdf_varget,cdfid,rhum_id,rhum
    if omega_id ge 0 then ncdf_varget,cdfid,omega_id,omega
;
;  Close the netcdf file
;
    ncdf_close,cdfid
;
;  break up the arrays according to domain
;  we only want to use domain 1
;
    press=transpose(press)  ;change to time,levels,domains
    result=where(press ne -99999, count)
    if count ne 0 then press[result]=press[result]/100.  ;convert from pa to mb
    press=reverse(press,2)  ;reverse rows
    press_1=press[*,*,0] 
    press_1=reform(press_1)

    cfrac=transpose(cfrac)  ;time,levels,domains
    cfrac=reverse(cfrac,2)
    cfrac_1=cfrac[*,*,0]
    cfrac_1=reform(cfrac_1)

    uwind=transpose(uwind)  ;time,levels,domains
    uwind=reverse(uwind,2)
    uwind_1=uwind[*,*,0]
    uwind_1=reform(uwind_1)

    vwind=transpose(vwind)  ;time,levels,domains
    vwind=reverse(vwind,2)
    vwind_1=vwind[*,*,0]
    vwind_1=reform(vwind_1)

    temp=transpose(temp)  ;time,levels,domains
    temp=reverse(temp,2)
    temp_1=temp[*,*,0]
    temp_1=reform(temp_1)

    shum=transpose(shum)  ;time,levels,domains
    shum=reverse(shum,2)
    shum_1=shum[*,*,0]
    shum_1=reform(shum_1)

    liq=transpose(liq)  ;time,levels,domains
    liq=reverse(liq,2)
    liq_1=liq[*,*,0]
    liq_1=reform(liq_1)

    ice=transpose(ice)  ;time,levels,domains
    ice=reverse(ice,2)
    ice_1=ice[*,*,0]
    ice_1=reform(ice_1)

    rhum=transpose(rhum)  ;time,levels,domains
    rhum=reverse(rhum,2)
    rhum_1=rhum[*,*,0]
    rhum_1=reform(rhum_1)

    omega=transpose(omega)  ;time,levels,domains
    omega=reverse(omega,2)
    omega_1=omega[*,*,0]
    omega_1=reform(omega_1)
;
;  If it is the first file
;
    if i eq 0 then begin  
      com_btto=long(base_time)+long(time_offset)
      com_press_1=press_1
      com_cfrac_1=cfrac_1
      com_uwind_1=uwind_1
      com_vwind_1=vwind_1
      com_temp_1=temp_1
      com_shum_1=shum_1
      com_liq_1=liq_1
      com_ice_1=ice_1
      com_rhum_1=rhum_1
      com_omega_1=omega_1
;
;  It if is not the first file
;
    endif else begin  
      com_btto=[long(com_btto),long(base_time)+long(time_offset)]
      com_press_1=[com_press_1,press_1]
      com_cfrac_1=[com_cfrac_1,cfrac_1]
      com_uwind_1=[com_uwind_1,uwind_1]
      com_vwind_1=[com_vwind_1,vwind_1]
      com_temp_1=[com_temp_1,temp_1]
      com_shum_1=[com_shum_1,shum_1]
      com_liq_1=[com_liq_1,liq_1]
      com_ice_1=[com_ice_1,ice_1]
      com_rhum_1=[com_rhum_1,rhum_1]
      com_omega_1=[com_omega_1,omega_1]
    endelse
  endif  ;end of not a gif file
endfor  ;end of loop through matching files
;popd  ;return to program directory
;
;  Redefine arrays
;
press_1=com_press_1
cfrac_1=com_cfrac_1
uwind_1=com_uwind_1
vwind_1=com_vwind_1
temp_1=com_temp_1
shum_1=com_shum_1
liq_1=com_liq_1
ice_1=com_ice_1
rhum_1=com_rhum_1
omega_1=com_omega_1
;
;  Get parameters of new arrays
;
base_time=long(com_btto[0])  ;base time
time_offset=long(com_btto)-long(base_time)
numtimes=n_elements(time_offset)
;
;  Calculate height from pressure
;
if month eq '03' then pressure=press_1[fix(numtimes/2),*] else  pressure=press_1[0,*]
To=288  ;Kelvins
Po=1013.25  ;mb
G=6.50  ;deg km-1
G=G/1000.  ;deg m-1
grav=9.81  ;m/s
R=287  ;J/deg/kg
exponent=R*G/grav
bracket=1-(pressure/Po)^exponent
height=To/G*bracket
height=height/1000.  ;convert to km
;
;  Create new regular height grid
;
if site eq 'nsa' then numheights=61 else numheights=100
newheight=findgen(numheights)*0.25
;
;  Interpolate 2D grids to newheights grid
;
newpress=fltarr(numtimes,numheights)
newcfrac=fltarr(numtimes,numheights)
newuwind=fltarr(numtimes,numheights)
newvwind=fltarr(numtimes,numheights)
newtemp=fltarr(numtimes,numheights)
newshum=fltarr(numtimes,numheights)
newliq=fltarr(numtimes,numheights)
newice=fltarr(numtimes,numheights)
newrhum=fltarr(numtimes,numheights)
newomega=fltarr(numtimes,numheights)
;
;  Fixes the march problem in this loop
;  The bottom levels of March for the first 10 days are
;  -99999.  This messes up the interpolation.  Thus 
;  the arrays have to be restricted to only good values
;  before interpolation.
;
for t=0,numtimes-1 do begin

  tempo=press_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newpress[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newpress[t,result]=-99999
  endif else newpress[t,*]=-99999

  tempo=cfrac_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newcfrac[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newcfrac[t,result]=-99999
  endif else newcfrac[t,*]=-99999

  tempo=uwind_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newuwind[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newuwind[t,result]=-99999
  endif else newuwind[t,*]=-99999

  tempo=vwind_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newvwind[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newvwind[t,result]=-99999
  endif else newvwind[t,*]=-99999

  tempo=temp_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newtemp[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newtemp[t,result]=-99999
  endif else newtemp[t,*]=-99999

  tempo=shum_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newshum[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newshum[t,result]=-99999
  endif else newshum[t,*]=-99999

  tempo=liq_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newliq[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newliq[t,result]=-99999
  endif else newliq[t,*]=-99999

  tempo=ice_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newice[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newice[t,result]=-99999
  endif else newice[t,*]=-99999

  tempo=rhum_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newrhum[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newrhum[t,result]=-99999
  endif else newrhum[t,*]=-99999

  tempo=omega_1[t,*]
  result=where(tempo ne -99999,count)
  if count ne 0 then begin
    tempo=tempo[result]
    tempheight=height[0:count-1]
    newomega[t,*]=interpol(tempo,tempheight,newheight)
    result=where(newheight gt tempheight[n_elements(tempheight)-1],count) 
    if count ne 0 then newomega[t,result]=-99999
  endif else newomega[t,*]=-99999
;	newpress[t,*]=interpol(press_1[t,*],height,newheight)
;	newcfrac[t,*]=interpol(cfrac_1[t,*],height,newheight)
;	newuwind[t,*]=interpol(uwind_1[t,*],height,newheight)
;	newvwind[t,*]=interpol(vwind_1[t,*],height,newheight)
;	newtemp[t,*]=interpol(temp_1[t,*],height,newheight)
;	newshum[t,*]=interpol(shum_1[t,*],height,newheight)
;	newliq[t,*]=interpol(liq_1[t,*],height,newheight)
;	newice[t,*]=interpol(ice_1[t,*],height,newheight)
;	newrhum[t,*]=interpol(rhum_1[t,*],height,newheight)
;	newomega[t,*]=interpol(omega_1[t,*],height,newheight)
endfor
;
;  Set negative values to 0 in data where 0's are nonsensical.
;  These negative values result from the interpolation
;  Set values > 1 to 0 in data where >1 is nonsensical.
;  These large values result from the interpolation.
;
result=where(newliq lt 0,count)
if count ne 0 then newliq[result]=0.0
result=where(newice lt 0,count)
if count ne 0 then newice[result]=0.0
result=where(newshum lt 0,count)
if count ne 0 then newshum[result]=0
result=where(newcfrac lt 0 or newcfrac gt 1,count)
if count ne 0 then newcfrac[result]=0.0
result=where(newrhum lt 0 or newrhum gt 1,count)
if count ne 0 then newrhum[result]=0.0
;
;  Zip up the data files and make sure file permissions
;  and ownerships are okay
;
unix_command='gzip -fq '+arm_file  
spawn,unix_command
unix_command='chmod 664 '+arm_file
spawn,unix_command
unix_command='chgrp met-mace '+arm_file
spawn,unix_command
;
;  Return to the original program directory
;
popd

return
end