PRO Read_reference_ERA, input_file, $
                        pressure, $
                        temperature, $
                        watervapor, $  ; specific humidity
                        vapor_pressure, $  ; for refractivity
                        Geopotential_height

ncfileid = NCDF_OPEN(input_file)
NCDF_VARGET, ncfileid, 'T_GDS4_ISBL', temperature ; K, [512, 256, 37], float
NCDF_VARGET, ncfileid, 'Q_GDS4_ISBL', specific ; kg/kg,[512, 256, 37], float

NCDF_VARGET, ncfileid, 'Z_GDS4_ISBL', Geopotential_height  ; m**2 s**-2, [512, 256, 37], float
NCDF_VARGET, ncfileid, 'g4_lat_1', latd_single  ; 256

;NCDF_VARGET, ncfileid, 'lv_ISBL0', pressure ; 37, hPa
NCDF_CLOSE, ncfileid

temperature = double(temperature)
specific = double(specific)    

;print,string(10B),input_file    
;print,'size(specific): ',size(specific)
;print,'size(pressure): ',size(pressure)                         

; convert kg/kg to g/kg
watervapor=specific*1000.

; convert specific to watervapor pressure
vapor_pressure=[specific]
for pressure_loop = 0, n_elements(pressure)-1 do begin
  vapor_pressure[*,*,pressure_loop] = pressure[pressure_loop]*specific[*,*,pressure_loop]/(0.622+0.378*specific[*,*,pressure_loop])
endfor

END


FUNCTION SUB_INPUT_ERA, pdatestr_dashed, $
                        ; -------
                        lon=lon_era, $
                        lat=lat_era, $  
                        pressure=pressure_era, $ 
                        ; -------
                        Geopotential_height_all, $
                        Temperature_all, $
                        watervapor_all, $
                        ref=refractivity_all, $
                        dryt=dryTemperature_all

  ;print,'pdatestr_dashed',pdatestr_dashed
  pdatestr_array=strsplit(pdatestr_dashed,'-',/extract)
  pyearstr=pdatestr_array[0]
  pmonthstr=pdatestr_array[1]
  pdaystr=pdatestr_array[2]
  pdatestr = pyearstr+pmonthstr+pdaystr

  caldat,julday(pmonthstr,fix(pdaystr)+1,pyearstr), pmonthstr_tomorrow, pdaystr_tomorrow, pyearstr_tomorrow
  pyearstr_tomorrow = string(pyearstr_tomorrow, format='(i4)')
  pdatestr_tomorrow = pyearstr_tomorrow+string(pmonthstr_tomorrow, format='(i2.2)')+string(pdaystr_tomorrow, format='(i2.2)')

  ; ----------------------------
  spawn, 'hostname -s', hostname

  if strmid(hostname, 0, 3) eq 'rhw' then begin
    folder= '/data/data504/xzhou/'

  endif else if hostname eq 'gnssro' then begin
    folder= '/data/xinjiaz/'

  endif else if hostname eq 'starro' then begin  
    folder= '/GNSSRO_data/xinjiaz/'

  endif else if hostname eq 'gpsmet' then begin  
    folder= '/GNSSRO_data/xinjiaz/'

  endif
  
  path = folder+'data_ERA/ERAi_rda/'

  L4_patten_today    = path+pyearstr+'/ei.oper.an.pl.regn128sc.'+pdatestr+'??.nc'
  L4_patten_tomorrow = path+pyearstr_tomorrow+'/ei.oper.an.pl.regn128sc.'+pdatestr_tomorrow+'00.nc'
  L4_patten = [L4_patten_today, L4_patten_tomorrow]
  L4_FileList = FILE_SEARCH(L4_patten, COUNT=L4_patten_Count)

  if L4_patten_Count ne 5 then begin
    print, 'there is no enough ERA files in '+pdatestr_dashed
    temperature_all=!VALUES.F_NAN
    watervapor_all=!VALUES.F_NAN
    Geopotential_height_all=!VALUES.F_NAN
    dryTemperature_all=!VALUES.F_NAN
    refractivity_all=!VALUES.F_NAN

    return, 0
  endif else begin
  
    restore,path+'dim.sav'
    ;lon_era, [ 0 ~ 360]
    ;lat_era, [90 ~ -90]
    ;pressure_era, [1 ~ 1000], 37 levels, hPa

    temperature_all= []
    watervapor_all = []
    vapor_pressure_all = []
    Geopotential_height_all=[]

    for file_loop = 0, L4_patten_Count-1 do begin

      ; --------------------
      print,L4_FileList[file_loop]
      Read_reference_ERA, L4_FileList[file_loop], pressure_era, temperature, watervapor, vapor_pressure, Geopotential_height

      temperature = [temperature, temperature[0, *, *]]  ; add one more column for longitude
      watervapor  = [watervapor, watervapor[0, *, *]]
      vapor_pressure  = [vapor_pressure, vapor_pressure[0, *, *]]
      Geopotential_height = [Geopotential_height, Geopotential_height[0, *, *]]
  
      temperature_all = [temperature_all, reform(temperature, [1, 513, 256, 37])]
      watervapor_all = [watervapor_all, reform(watervapor, [1, 513, 256, 37])]
      vapor_pressure_all = [vapor_pressure_all, reform(vapor_pressure, [1, 513, 256, 37])]
      Geopotential_height_all= [Geopotential_height_all, reform(Geopotential_height, [1, 513, 256, 37])]

    endfor

    ; for interpolation, mv 'vertical' to left most, because interpolation starting from rightmost
    Temperature_all = transpose(temperature_all, [3, 0, 1, 2])  ; 37:5:513:256, height:time:lon:lat
    watervapor_all  = transpose(watervapor_all, [3, 0, 1, 2])
    vapor_pressure_all  = transpose(vapor_pressure_all, [3, 0, 1, 2])
    Geopotential_height_all = transpose(Geopotential_height_all, [3, 0, 1, 2])

    ; rebuild pressure for refractivity and dry temperature
    size_dim = size(Temperature_all, /dim)
    pressure_all = dblarr(size_dim)
    for ii = 0, size_dim[0]-1 do pressure_all[ii,*,*,*]=pressure_era[ii]

    ; calculate refractivity
    if ARG_PRESENT(refractivity_all) then refractivity_all = 77.6D*pressure_all/Temperature_all + 3.73e5 * vapor_pressure_all / Temperature_all / Temperature_all

    ; calculate dry temperature
    if ARG_PRESENT(dryTemperature_all) then dryTemperature_all = 77.6D*pressure_all/(77.6D*pressure_all/Temperature_all + 3.73e5 * vapor_pressure_all / Temperature_all / Temperature_all)

    return, 2
  endelse
END