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

ncfileid = NCDF_OPEN(input_file)
NCDF_VARGET, ncfileid, 't', temperature ; K, [1440, 721, 37, 4], float
NCDF_VARGET, ncfileid, 'q', specific   ; kg/kg, [1440, 721, 37, 4], float
NCDF_VARGET, ncfileid, 'z', Geopotential_height  ; m**2 s**-2, [1440, 721, 37, 4], float
NCDF_CLOSE, ncfileid

; [1000 ~ 1] --> [1 ~ 1000]
temperature = reverse(temperature,3)
specific = reverse(specific,3)
Geopotential_height = reverse(Geopotential_height,3)

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

if keyword_set(for_tomorrow) then begin ; only keep the first time
  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_cds_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_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 'gpsmet' then begin  
    folder= '/ROSDC_data3/xinjiaz/'

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

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

  endif else if hostname eq 'rosdc' then begin  
    folder= '/data3/xinjiaz/'

  endif

  path = folder+'data_ERA/ERA5_cds_v2/'

  ;                                                           e5.oper.an.pl.128_20240411.nc
  L4_file_today    = path+pyearstr         +'/e5.oper.an.pl.128_'+pdatestr         +'.nc'
  L4_file_tomorrow = path+pyearstr_tomorrow+'/e5.oper.an.pl.128_'+pdatestr_tomorrow+'.nc'

  print, ''
  if ~file_test(L4_file_today) or ~file_test(L4_file_tomorrow) then begin
    print, 'there is no enough ERA5_cds_v2 files in '+pdatestr_dashed
    if ~file_test(L4_file_today)    then print, L4_file_today
    if ~file_test(L4_file_tomorrow) then print, L4_file_tomorrow

    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'
    ;height_era5, 37,  [47 ~ 0 ] km
    ;lon_era5, 1440+1, [0 ~ 359.75]+360.00
    ;lat_era5, 721, [ 90 ~ -90]
    ;pressure_era5, 37,  [1 ~ 1000] hPa (v1) <-- ICVS only accept high to surface. Reverse variables in "Read_reference_ERA5_cds_v2"
    ;pressure_era5, 37,  [1000 ~ 1] hPa (v2)

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

    ; --------------------
    print, L4_file_today
    Read_reference_ERA5_cds_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, 4] -> [1441, 721, 37, 4]
    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,4] -> [4, 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_cds_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]])

    ; --------------------
    temperature_all= [Temperature, temperature_tomorrow]  ; combine today and tomorrow
    specific_all = [specific, specific_tomorrow]   ; [4, 1441, 721, 37] -> [5, 1441, 721, 37]
    vapor_pressure_all = [vapor_pressure, vapor_pressure_tomorrow]   
    Geopotential_height_all= [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])   ; [5, 1441, 721, 37] -> [37, 5, 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