PRO Read_reference_ERA5_ads_v2,   input_file, $
                          pressure, $
                          temperature, $
                          specific, $  ; specific humidity
                          vapor_pressure, $  ; for refractivity
                          Geopotential_height, $
                          for_yesterday=for_yesterday, $
                          for_tomorrow=for_tomorrow

ncfileid = NCDF_OPEN(input_file)
NCDF_VARGET, ncfileid, 't', temperature ; K, [1440, 721, 37, 2], float
NCDF_attget, ncfileid, 't', 'scale_factor', scale_factor
NCDF_attget, ncfileid, 't', 'add_offset',  add_offset
temperature = temperature*scale_factor+add_offset

NCDF_VARGET, ncfileid, 'q', specific   ; kg/kg, [1440, 721, 37, 2], float
NCDF_attget, ncfileid, 'q', 'scale_factor', scale_factor
NCDF_attget, ncfileid, 'q', 'add_offset',  add_offset
specific = specific*scale_factor+add_offset

NCDF_VARGET, ncfileid, 'z', Geopotential_height  ; m**2 s**-2, [1440, 721, 37, 2], float
NCDF_attget, ncfileid, 'z', 'scale_factor', scale_factor
NCDF_attget, ncfileid, 'z', 'add_offset',  add_offset
Geopotential_height = Geopotential_height*scale_factor+add_offset
NCDF_CLOSE, ncfileid

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

if keyword_set(for_yesterday) then begin ; only keep the last time (18:00)
  temperature = reform(temperature[*, *, *, 1])
  specific = reform(specific[*, *, *, 1])
  Geopotential_height = reform(Geopotential_height[*, *, *, 1])
endif else if keyword_set(for_tomorrow) then begin ; only keep the first time (06:00)
  temperature = reform(temperature[*, *, *, 0])
  specific = reform(specific[*, *, *, 0])
  Geopotential_height = reform(Geopotential_height[*, *, *, 0])
endif

; -------------------
; convert m^2/s^2 to gpm
g_wmo = 9.80665D     ; reference gravity
Geopotential_height /= g_wmo

; -------------------
; calculate watervapor pressure, for refractivity
vapor_pressure=dblarr(size(specific, /dim))
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

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

END


FUNCTION SUB_INPUT_ERA5_ads_v2,   pdatestr_dashed, $
                          ; -------
                          lon=lon_era5, $
                          lat=lat_era5, $  
                          pressure=pressure_era5, $ 
                          ; -------
                          Geopotential_height_all, $
                          Temperature_all, $
                          specific_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_yesterday, pdaystr_yesterday, pyearstr_yesterday
  pyearstr_yesterday = string(pyearstr_yesterday, format='(i4)')
  pmonthstr_yesterday = string(pmonthstr_yesterday, format='(i2.2)')
  pdatestr_yesterday = pyearstr_yesterday+pmonthstr_yesterday+string(pdaystr_yesterday, format='(i2.2)')

  caldat,julday(pmonthstr,fix(pdaystr)+1,pyearstr), pmonthstr_tomorrow, pdaystr_tomorrow, pyearstr_tomorrow
  pyearstr_tomorrow = string(pyearstr_tomorrow, format='(i4)')
  pmonthstr_tomorrow = string(pmonthstr_tomorrow, format='(i2.2)')
  pdatestr_tomorrow = pyearstr_tomorrow+pmonthstr_tomorrow+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/ERA5_forecast/'

  ; 2022/e5.oper.fc.pl.128_20221001.nc
  L4_file_yesterday = file_search(path+pyearstr_yesterday+'/e5.oper.fc.pl.128_'+pdatestr_yesterday+'.nc', count=count_yesterday)
  L4_file_today     = file_search(path+pyearstr          +'/e5.oper.fc.pl.128_'+pdatestr          +'.nc', count=count_today)
  L4_file_tomorrow  = file_search(path+pyearstr_tomorrow +'/e5.oper.fc.pl.128_'+pdatestr_tomorrow +'.nc', count=count_tomorrow)

  print, ''
  if count_today eq 0 or count_yesterday eq 0 or count_tomorrow eq 0 then begin
    print, 'there is no enough ERA5_forecast files in '+pdatestr_dashed
    print, path+pyearstr_yesterday+'/e5.oper.fc.pl.128_'+pdatestr_yesterday+'.nc'
    print, path+pyearstr          +'/e5.oper.fc.pl.128_'+pdatestr          +'.nc'
    print, path+pyearstr_tomorrow +'/e5.oper.fc.pl.128_'+pdatestr_tomorrow +'.nc'

    temperature_all=!VALUES.F_NAN
    specific_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_era5, 1440, [0 ~ 359.75]
    ;lat_era5, 721, [ 90 ~ -90]
    ;pressure_era5, 37,  [1 ~ 1000 ]

    temperature_all= []
    specific_all = []
    vapor_pressure_all = []
    Geopotential_height_all=[]

    ; --------------------
    print, L4_file_yesterday
    Read_reference_ERA5_ads_v2, L4_file_yesterday, $
                          pressure_era5, $
                          temperature_yesterday, $
                          specific_yesterday, $
                          vapor_pressure_yesterday, $
                          Geopotential_height_yesterday, $
                          for_yesterday=1

    temperature_yesterday = [temperature_yesterday, temperature_yesterday[0, *, *]]  ; add one more column for longitude
    specific_yesterday  = [specific_yesterday, specific_yesterday[0, *, *]]          ; [1440, 721, 37] -> [1441, 721, 37]
    vapor_pressure_yesterday  = [vapor_pressure_yesterday, vapor_pressure_yesterday[0, *, *]] 
    Geopotential_height_yesterday = [Geopotential_height_yesterday, Geopotential_height_yesterday[0, *, *]]
  
    size_dim = size(temperature_yesterday, /dim)
    temperature_yesterday = reform(temperature_yesterday, [1, size_dim[0], size_dim[1], size_dim[2]])  ; add one column, for combining data
    specific_yesterday = reform(specific_yesterday, [1, size_dim[0], size_dim[1], size_dim[2]])
    vapor_pressure_yesterday = reform(vapor_pressure_yesterday, [1, size_dim[0], size_dim[1], size_dim[2]])
    Geopotential_height_yesterday = reform(Geopotential_height_yesterday, [1, size_dim[0], size_dim[1], size_dim[2]])
    ; --------------------
    print, L4_file_today
    Read_reference_ERA5_ads_v2, L4_file_today, $
                          pressure_era5, $
                          temperature, $
                          specific, $
                          vapor_pressure, $
                          Geopotential_height

    temperature = [temperature, temperature[0, *, *, *]] ; add one more column for longitude
    specific  = [specific, specific[0, *, *, *]]         ; [1440, 721, 37, 2] -> [1441, 721, 37, 2]
    vapor_pressure  = [vapor_pressure, vapor_pressure[0, *, *, *]]   
    Geopotential_height = [Geopotential_height, Geopotential_height[0, *, *, *]]
  
    Temperature = transpose(temperature, [3, 0, 1, 2])  ; mv time to first column, for combining with data in tomorrow
    specific  = transpose(specific, [3, 0, 1, 2])       ; [1441, 721, 37,2] -> [2, 1441, 721, 37]
    vapor_pressure  = transpose(vapor_pressure, [3, 0, 1, 2])   
    Geopotential_height = transpose(Geopotential_height, [3, 0, 1, 2])

    ; --------------------
    print, L4_file_tomorrow
    Read_reference_ERA5_ads_v2, L4_file_tomorrow, $
                          pressure_era5, $
                          temperature_tomorrow, $
                          specific_tomorrow, $
                          vapor_pressure_tomorrow, $
                          Geopotential_height_tomorrow, $
                          for_tomorrow=1

    temperature_tomorrow = [temperature_tomorrow, temperature_tomorrow[0, *, *]]  ; add one more column for longitude
    specific_tomorrow  = [specific_tomorrow, specific_tomorrow[0, *, *]]          ; [1440, 721, 37] -> [1441, 721, 37]
    vapor_pressure_tomorrow  = [vapor_pressure_tomorrow, vapor_pressure_tomorrow[0, *, *]] 
    Geopotential_height_tomorrow = [Geopotential_height_tomorrow, Geopotential_height_tomorrow[0, *, *]]
  
    size_dim = size(temperature_tomorrow, /dim)
    temperature_tomorrow = reform(temperature_tomorrow, [1, size_dim[0], size_dim[1], size_dim[2]])  ; add one column, for combining data
    specific_tomorrow = reform(specific_tomorrow, [1, size_dim[0], size_dim[1], size_dim[2]])
    vapor_pressure_tomorrow = reform(vapor_pressure_tomorrow, [1, size_dim[0], size_dim[1], size_dim[2]])
    Geopotential_height_tomorrow = reform(Geopotential_height_tomorrow, [1, size_dim[0], size_dim[1], size_dim[2]])

    ; --------------------
    ; combine yesterday, today and tomorrow
    ; [2, 1441, 721, 37] -> [4, 1441, 721, 37]
    temperature_all= [Temperature_yesterday, Temperature, temperature_tomorrow] 
    specific_all = [specific_yesterday, specific, specific_tomorrow]   
    vapor_pressure_all = [vapor_pressure_yesterday, vapor_pressure, vapor_pressure_tomorrow]   
    Geopotential_height_all= [Geopotential_height_yesterday, Geopotential_height, Geopotential_height_tomorrow]  

    ; --------------------
    Temperature_all = transpose(temperature_all, [3, 0, 1, 2])  ; mv 'vertical' to left most, because interpolation starting from rightmost
    specific_all  = transpose(specific_all, [3, 0, 1, 2])   ; [3, 1441, 721, 37] -> [37, 3, 1441, 721]
    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_era5[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.6*pressure_all/(77.6D*pressure_all/Temperature_all + 3.73e5 * vapor_pressure_all / Temperature_all / Temperature_all)

    return, 2
  endelse
END