PRO abi_nacp_sipnet,prefix,indir=indir,outdir=outdir
;Read "prefix"forcing.nc
;convert to Sipnet "prefix".clim and .spd file

;input: NACP met file

;output: SIPNET ASCII files

;sipnet.clim
;Location Year Day Hour Timestep(fday)
;Tair(C) Tsoil(C) PAR(ein/m2/timestep) Precip(mm/timestep)
;VPD (PA) VPDS (Pa) Vpress (Pa) U (m/s) S_moist(%sat)

;sipnet.spd
;Year_start jd_start steps per day #number of days -1

;Set up file names and default folders 
  IF n_elements(indir) EQ 0 THEN indir = expand_path('/Users/adesai/Desktop')+'/'
  IF n_elementS(outdir) EQ 0 THEN outdir = expand_path('/Users/adesai/Desktop')+'/'
  IF n_elements(prefix) EQ 0 THEN prefix = 'US-WCr'
  infile = indir + prefix + 'forcing.nc'
  climfile = outdir + prefix + '.clim'
  spdfile = outdir + prefix + '.spd'

;Test for file read/write
  IF file_test(infile,/read) AND file_test(outdir,/write) THEN BEGIN
    nc = ncdf_open(infile)
    
;From the attributes of the time variable, find the start year and day
    ncdf_attget,nc,'t','time_origin',time_origin
    time_origin = strmid(string(time_origin),0,10)
    year = strmid(time_origin,0,4)
    month = strmid(time_origin,5,2)
    day = strmid(time_origin,8,2)
    DOY = dy_to_jd(year+month+day)
    
;From attributes of timestep and dimensions, find time step and number of days
    ncdf_attget,nc,'timestp','tstep_sec',tstep_sec
    steps_per_day = 86400l/long(tstep_sec)
    ncdf_diminq,nc,0,tdimname,tdim ;first dimension is time
    ndays = tdim / steps_per_day
    
;Build string to write to spd file
    spd = year + ' ' + string(long(DOY),format='(i0)') + ' ' + string(steps_per_day,format='(i0)') + ' #' + string(ndays,format='(i0)') + ' -1'

;Build time output for the clim file (year day hour timestep)

;Create a Julian time array, convert to calendar time, convert to
;string, and then to day of year
;Convoluted but very fast and general - handles leap years et al gracefully

    times_j = TIMEGEN(tdim, UNITS='Seconds', STEP_SIZE=tstep_sec, START=JULDAY(long(month),long(day),long(year)))
    caldat,times_j,mon,day,yr,hr,mn,sec
    time_str = string(yr,format='(i4.4)')+string(mon,format='(i2.2)')+string(day,format='(i2.2)')
    jd = dy_to_jd(time_str)
    hour = float(hr)+(float(mn)/60.)
    timestep = float(tstep_sec) / 3600.
    
;Build met variables

;Air temp - convert to C
    ncdf_varget,nc,'Tair',Tair
    Tair -= 273.15

;Soil temp - model from air temperature
    Tsoil = smooth(tair,steps_per_day*7.0,/edge) > 0

;PAR - model from SW
    ncdf_varget,nc,'SWdown',SW
    PAR = SW*2.21  ;4.43 umol/J * 0.5
    PAR = PAR * float(tstep_sec) / 1e6 ;convert to einsteins/m2/timestep

;Precip(mm/timestep)
    ncdf_varget,nc,'Rainf',precip  ;kg/m2/s = mm/s
    precip = precip * float(tstep_sec)  ;mm/timestep
    
;VPD (Pa),VPDS (Pa), Vpress (Pa)
    ncdf_varget,nc,'Qair',Qair ;specific humidity kg/kg
    ncdf_varget,nc,'Psurf',Press ;pressure Pa
    Es = 100. * 6.112 * exp( (17.67 * Tair) / (243.5 + Tair)) ;saturation vapor pressure
    Es_soil = 100. * 6.112 * exp( (17.67 * Tsoil) / (243.5 + Tsoil))
    Ea = (Qair * Press)/(0.622 + (1 - 0.622)* Qair) ;vapor pressure
    vpd = es-ea > 0 ;VPD, correct for a few stray negative VPD
    vpd_s = es_soil - ea ;relative soil vapor deficit assuming soil pore space is saturated

;Wind speed (m/s)
    ncdf_varget,nc,'Wind',Wind

;S_moist(%sat) - set to fixed value 
    Smoist = 1.0
    
;Write the files

;sipnet.clim
;Location Year Day Hour Timestep(fday)
;Tair(C) Tsoil(C) PAR(ein/m2/timestep) Precip(mm/timestep)
;VPD (PA) VPDS (Pa) Vpress (Pa) U (m/s) S_moist(%sat)

;sipnet.spd
;Year_start jd_start steps per day #number of days -1

    openw,fl,climfile,/get_lun
    FOR i = 0,n_elements(yr)-1 DO BEGIN 
      printf,fl,0,yr[i],jd[i],hour[i],timestep,Tair[i],Tsoil[i],$
             PAR[i],Precip[i],VPD[i],VPD_S[i],Ea[i],Wind[i],Smoist,$
             format='(i1," ",i4," ",i3," ",f6.1," ",f15.8,9f15.3)'
    ENDFOR 
    free_lun,fl
    
    openw,fl,spdfile,/get_lun
    printf,fl,spd,format='(a0)'
    free_lun,fl
  ENDIF ELSE print,'File read/write error'
  
  stop
  
END