PRO readnarr_Acme,yr=yr

;use closest function to map lat
;for ACME 37 to 42 =   37.125 to 42 (15 pts)
;for ACME 111 to 105 = -110.875 to -104.875 (17 pts)
;all time

;vars
;dpt2m (2m dew point)
;rh2m
;tmp2m (k)
;tsoil0_10cm (k)
;soill0_10cm
;ugrd10m vgrd10m
;apcp (kg/m2)
;hpbl (m)
;presnsfc (sfc pressure)
;dswrfsfc (w/m2 shortwave)

  vars = ['tmp2m','dpt2m','rh2m','tsoil0_10cm','soill0_10cm','ugrd10m','vgrd10m','apcp','hpbl','presnsfc','dswrfsfc','snod','snowc','albdosfc','veg']

;also try to get 700mb vertical velocity (maybe not 700 in mountain,
;near 300 above ground?)
;actually vvelprs ugrdprs vgrdprs 

;snod (snow depth), snowc (snow cover), albdosfc, veg, 
;high/low/mid clouds

;1/1/07-12/31/07
;AWIP32-fixed
;vgtypsfc
;hgtsfc

;for each day
;read the grids 

;vars is an array of vars - copy from structs to array

  IF n_elements(yr) EQ 0 THEN yr = 2007

  lats = numgen(37.125,42.0,numvals=14)
  lons = numgen(-110.875,-104.875,numvals=17)
  tms_utc = findgen(8)*3
  airt = fltarr(366,17,14,8)
  dewpt = airt
  rh = airt
  soilt = airt
  soilq = airt
  u = airt
  v = airt
  precip = airt
  zi = airt
  press = airt
  sw = airt
  snowd = airt
  snowc = airt
  albedo = airt
  vegcov = airt

;read fixed fields - collapse to original grid
  print,'Reading fixed data'
  vv = ['hgtsfc','vgtypsfc']
  fx = readnarr(2007,1,1,vv,lats=[37,42],lons=[-111,-104.875],/fixed)
  elevhi = fx.hgtsfc.data 
  elev = rebin(elevhi,17,14)
  veghi = fx.vgtypsfc.data
  veg = rebin(veghi,17,14,/sample)
  
;1 = EBF, 2 = DBF, 3 = Mixed, 4 = ENF, 5 = DNF, 6 = C4, 7 = Grass, 8 =
;shrub, NOT IN ACME: 9 = bare soil, 10 = tundra, 11 = desert, 12 = ag

;for modeling, 
;1,2,3,4,5 = Niwot, 6,7,8 GLEES or nothing  (< 5.5

;timezone = UTC-7 MST
;0 UTC = 5pm previous day

;read dailys
  FOR mo = 1,12 DO BEGIN
    FOR day = 1,days_in_month(mo,y=yr) DO BEGIN 
      doy = dy_to_jd(string(yr,format='(i4.4)')+string(mo,format='(i2.2)')+string(day,format='(i2.2)'))
      print,'On Doy',doy
      t = readnarr(yr,mo,day,vars,lats=[37,42],lons=[-111,-105])

; vars = ['tmp2m','dpt2m','rh2m','tsoil0_10cm','soill0_10cm','ugrd10m','vgrd10m','apcp','hpbl','presnsfc','dswrfsfc']
      airt[doy-1,*,*,*] = t.tmp2m.data
      dewpt[doy-1,*,*,*] = t.dpt2m.data
      rh[doy-1,*,*,*] = t.rh2m.data
      soilt[doy-1,*,*,*] = t.tsoil0_10cm.data
      soilq[doy-1,*,*,*] = t.soill0_10cm.data
      u[doy-1,*,*,*] = t.ugrd10m.data
      v[doy-1,*,*,*] = t.vgrd10m.data
      precip[doy-1,*,*,*] = t.apcp.data
      zi[doy-1,*,*,*] = t.hpbl.data
      press[doy-1,*,*,*] = t.presnsfc.data
      sw[doy-1,*,*,*] = t.dswrfsfc.data
      snowd[doy-1,*,*,*] = t.snod.data
      snowc[doy-1,*,*,*] = t.snowc.data
      albedo[doy-1,*,*,*] = t.albdosfc.data
      vegcov[doy-1,*,*,*] = t.veg.data

    ENDFOR
    print,'Saving in month ',mo
    save,elevhi,elev,veghi,veg,airt,dewpt,rh,soilt,soilq,u,v,precip,zi,press,sw,snowd,snowc,albedo,vegcov,lats,lons,tms_utc,filename=mydocs()+'acme07/data/narr/acmenarr'+string(yr,format='(i4.4)')+'.temp.sav'
  ENDFOR

;read Jan 1 of next year
  t = readnarr(yr+1,1,1,vars,lats=[37,42],lons=[-111,-105])
  airt[doy,*,*,*] = t.tmp2m.data
  dewpt[doy,*,*,*] = t.dpt2m.data
  rh[doy,*,*,*] = t.rh2m.data
  soilt[doy,*,*,*] = t.tsoil0_10cm.data
  soilq[doy,*,*,*] = t.soill0_10cm.data
  u[doy,*,*,*] = t.ugrd10m.data
  v[doy,*,*,*] = t.vgrd10m.data
  precip[doy,*,*,*] = t.apcp.data
  zi[doy,*,*,*] = t.hpbl.data
  press[doy,*,*,*] = t.presnsfc.data
  sw[doy,*,*,*] = t.dswrfsfc.data
  snowd[doy,*,*,*] = t.snod.data
  snowc[doy,*,*,*] = t.snowc.data
  albedo[doy,*,*,*] = t.albdosfc.data
  vegcov[doy,*,*,*] = t.veg.data
;,'snod','snoc','albdosfc','veg'


;transpose vars
  airt = reform(transpose(airt,[3,0,1,2]),2928,17,14)
  dewpt = reform(transpose(dewpt,[3,0,1,2]),2928,17,14)
  rh = reform(transpose(rh,[3,0,1,2]),2928,17,14)
  soilt = reform(transpose(soilt,[3,0,1,2]),2928,17,14)
  soilq = reform(transpose(soilq,[3,0,1,2]),2928,17,14)
  u = reform(transpose(u,[3,0,1,2]),2928,17,14)
  v = reform(transpose(v,[3,0,1,2]),2928,17,14)
  precip = reform(transpose(precip,[3,0,1,2]),2928,17,14)
  zi = reform(transpose(zi,[3,0,1,2]),2928,17,14)
  press = reform(transpose(press,[3,0,1,2]),2928,17,14)
  sw = reform(transpose(sw,[3,0,1,2]),2928,17,14)
  albedo = reform(transpose(albedo,[3,0,1,2]),2928,17,14)
  snowc = reform(transpose(snowc,[3,0,1,2]),2928,17,14)
  snowd = reform(transpose(snowd,[3,0,1,2]),2928,17,14)
  vegcov = reform(transpose(vegcov,[3,0,1,2]),2928,17,14)

;timezone shift
  airt = airt[2:2921,*,*]
  dewpt = dewpt[2:2921,*,*]
  rh = rh[2:2921,*,*]
  soilt = soilt[2:2921,*,*]
  soilq = soilq[2:2921,*,*]
  u = u[2:2921,*,*]
  v = v[2:2921,*,*]
  precip = precip[2:2921,*,*]
  zi = zi[2:2921,*,*]
  press = press[2:2921,*,*]
  sw = sw[2:2921,*,*]
  snowd = snowd[2:2921,*,*]
  snowc = snowc[2:2921,*,*]
  vegcov = vegcov[2:2921,*,*]
  albedo = albedo[2:2921,*,*]

  tms = 1.0+(findgen(365.*8.)/8.)

;badval clean up - remove > 9e20 SKIP

;save vars
  save,elevhi,elev,veghi,veg,airt,dewpt,rh,soilt,soilq,u,v,precip,zi,press,sw,snowd,snowc,albedo,vegcov,lats,lons,tms,filename=mydocs()+'acme07/data/narr/acmenarr'+string(yr,format='(i4.4)')+'.sav'


  stop

;---
;sipnet conversion
;textfile for each gridpoint - label by lat/lon
;convert units - write file

END

PRO acme_convprism
;cellsize: .00833333333304444
;North Coordinate: 41.3666666
;South Coordinate: 36.0000000
;West Coordinate: -111.9666666
;East Coordinate: -104.8666666
;# Rows: 639
;# Columns: 853

  dr = mydocs()+'acme07/data/narr/prismcsv/'

;coords
  lats_acme = numgen(36.0d,41.36666666d,numvals=639)
  lons_acme = numgen(-111.9666666d,-104.86666666d,numvals=853)

;dem
  print,'DEM'
  dem_acme = (read_ascii(dr+'800mDEMEXCEL.csv',delim=',')).(0)
 
;tmax
  tmax_acme = fltarr(12,853,639)
  FOR i = 1,12 DO BEGIN
    print,'Tmax month ',i
    dta = (read_ascii(dr+'TMax2007ACMEExtract-'+string(i,format='(i0)')+'.csv',delim=',')).(0)
    tmax_acme[i-1,*,*] = dta
  ENDFOR 

  save,lats_acme,lons_acme,dem_acme,tmax_acme,filename=mydocs()+'acme07/data/narr/acmeprism_2007.sav'

;tmin
  tmin_acme = fltarr(12,853,639)
  FOR i = 1,12 DO BEGIN
    print,'Tmin month ',i
    dta = (read_ascii(dr+'TMin2007ACMEExtract-'+string(i,format='(i0)')+'.csv',delim=',')).(0)
    tmin_acme[i-1,*,*] = dta
  ENDFOR 

  save,lats_acme,lons_acme,dem_acme,tmax_acme,tmin_acme,filename=mydocs()+'acme07/data/narr/acmeprism_2007.sav'

END


;5 param files
;Niwot
;Deciduous
;C4/Shrub Grassland Shrub

;daily - scale
;use location for each 

;DAVE MONDAY: 5 param files

;stat downscale
;wilkes SDSM 

;priority: AGU spatial run


;.clim
;Location
;Year
;Day
;Hour
;Timestep (fraction day) 30 min = 0.02083333333 1 hr = 0.041666666667
;Tair - C - filledmet
;Tsoil - C - filledmet
;PAR - total einsteins/m2  = par*1800/1e6 - filledmet
;Precip mm = keep - filledmet?
;VPD Pa
;VPDsoil (Ts) - Pa = compute  (mr = mixing_ratio(tempC,pressMB,rh0-1,vpd=vpd kPA,Es=Es,Ea=Ea(mb) )
;Vpress (e) - Pa (1 mb = 100 Pa so Ea*100  1mb = .1Kpa so *0.1)
;Windspeed - m/s
;Soil wetness - fraction sat.

PRO acme_narrsipnet
;restore the NARR data
  restore,mydocs()+'acme07/data/narr/acmenarr2007.sav'
  restore,mydocs()+'acme07/data/narr/acmenarr_soilt2007.sav'

;reform and transpose
  airt = transpose(reform(airt,2920,17.*14.))-273.15
  soilt = transpose(reform(soilt25,2920,17.*14.))-273.15
  sw = transpose(reform(sw,2920,17.*14.))
  precip = transpose(reform(precip,2920,17.*14.))
  dewpt = transpose(reform(dewpt,2920,17.*14.))-273.15
  soilq = transpose(reform(soilq,2920,17.*14.))
  u = transpose(reform(u,2920,17.*14.))
  v = transpose(reform(v,2920,17.*14.))

;compute wind speed (u^2 + v^2 sqrt)
  ws = sqrt(u^2 + v^2)
;compute vp from Td
;compute vpd from Ta and vp
;compute vpd_S from Ts and vp
  vp = 100. * 6.112 * exp ( (17.67 * dewpt) / (243.5 + dewpt) )
  es = 100. * 6.112 * exp ( (17.67 * airt) / (243.5 + airt) )
  ess = 100. * 6.112 * exp ( (17.67 * soilt) / (243.5 + soilt) )
  vpd = (es-vp)>0
  vpds = ess-vp

;average_cols to daily (8)
;(sum for precip, avg for all others)

  airt_d = reform(transpose(average_cols(airt,8)),365.*17.*14.)
  soilt_d = reform(transpose(average_cols(soilt,8)),365.*17.*14.)
  sw_d = reform(transpose(average_cols(sw,8)),365.*17.*14.)
  precip_d = reform(transpose(average_cols(precip,8,/tot)),365.*17.*14.)
  soilq_d = reform(transpose(average_cols(soilq,8)),365.*17.*14.)
  ws_d = reform(transpose(average_cols(ws,8)),365.*17.*14.)
  vp_d = reform(transpose(average_cols(vp,8)),365.*17.*14.)
  vpd_d = reform(transpose(average_cols(vpd,8)),365.*17.*14.)
  vpds_d = reform(transpose(average_cols(vpds,8)),365.*17.*14.)

;Notes on units
;airt -> Tair  (conv K to C)
;soilt -> Tsoil (conv K to C)
;SW -> PAR (J/m2/s to mol/m2/day)  (~4)
;Precip -> Pre (kg/m2 * 1m3/1000kg * 1000mm/1m)  (no conv)
;VPD VPDs VP -> VPD VPDs VP in Pa
;ws -> wind in m/s
;sw-> soilwet in fraction (/100)

;location step from 0-237
  SW_TO_PAR = 2.03 * 86400./1e6
  clim = fltarr(14,365.*17.*14.)
  clim[0,*] = fix(findgen(365.*17.*14.)/365.)
  clim[1,*] = 2007.
  clim[2,*] = replicate_arr(findgen(365.)+1.,17.*14.)
  clim[3,*] = 0.
  clim[4,*] = 1.
  clim[5,*] = airt_d
  clim[6,*] = soilt_d
  clim[7,*] = sw_d * SW_TO_PAR 
  clim[8,*] = precip_d
  clim[9,*] = vpd_d
  clim[10,*] = vpds_d
  clim[11,*] = vp_d
  clim[12,*] = ws_d
  clim[13,*] = soilq_d

  stop

  write_xdf,mydocs()+'acme07/data/narr/acmeclim.xdf',clim
  form = '(i3," ",i4," ",i3," ",i1," ",i1," ",f14.3,f14.3,f14.3,f14.3,f14.3,f14.3,f14.3,f14.3,f14.3)'
  openw,fl,mydocs()+'acme07/data/narr/acme.clim',/get_lun
  printf,fl,clim,format=form
  free_lun,fl

  stop

;second file - veg
;Location Veg

;3 = 22  Decid and evergreen trees
;4 = 56  Evergreen needleleaf
;5 = 14  Deciduous needleleaf
;6 = 17  groundcover with tree/shrub
;7 = 17  groundcover only
;8 = 112 shrubs with perennial ground cover

;1J/2.02 mol
;J/s * 4.6 mol/J = mol/m2/s

;4.6 umol/m2/s / W/m2
;* 86400/1e6

END

;redo soil_temp at a deeper level
;tsoil10_40 cm


PRO readnarr_Acme_soilt,yr=yr
  vars = ['tsoil10_40cm']
  IF n_elements(yr) EQ 0 THEN yr = 2007
  lats = numgen(37.125,42.0,numvals=14)
  lons = numgen(-110.875,-104.875,numvals=17)
  tms_utc = findgen(8)*3
  soilt25 = fltarr(366,17,14,8)
  FOR mo = 1,12 DO BEGIN
    FOR day = 1,days_in_month(mo,y=yr) DO BEGIN 
      doy = dy_to_jd(string(yr,format='(i4.4)')+string(mo,format='(i2.2)')+string(day,format='(i2.2)'))
      print,'On Doy',doy
      t = readnarr(yr,mo,day,vars,lats=[37,42],lons=[-111,-105])
      soilt25[doy-1,*,*,*] = t.tsoil10_40cm.data
    ENDFOR
  ENDFOR

;read Jan 1 of next year
  t = readnarr(yr+1,1,1,vars,lats=[37,42],lons=[-111,-105])
  soilt25[doy,*,*,*] = t.tsoil10_40cm.data
  soilt25 = reform(transpose(soilt25,[3,0,1,2]),2928,17,14)
  soilt25 = soilt[2:2921,*,*]
;save vars
  save,soilt25,filename=mydocs()+'acme07/data/narr/acmenarr_soilt'+string(yr,format='(i4.4)')+'.sav'
  stop
END


PRO readnarr_point,lat,lon,yr=yr
;runs this from August 2005-October 2008

  vars = ['tmp2m','dpt2m','rh2m','tsoil10_40cm','soill0_10cm','ugrd10m','vgrd10m','apcp','hpbl','presnsfc','dswrfsfc','snod','snowc']

  IF n_elements(lat) EQ 0 THEN lat = 39.91
  IF n_elements(lon) EQ 0 THEN lon = -105.88
  IF n_elements(yr) EQ 0 THEN yr = 2007
  yrstr = string(yr,format='(i4.4)')
  tms_utc = findgen(8)*3
  airt = make_array(366,1,1,8,/float,value=nan())
  dewpt = airt
  rh = airt
  soilt = airt
  soilq = airt
  u = airt
  v = airt
  precip = airt
  zi = airt
  press = airt
  sw = airt
  snowd = airt
  snowc = airt

  endmo = 12
  startmo = 1

;if yr eq 2003 then startmo = 12
;if yr eq 2003 then restore,mydocs()+'acme07/data/narr/fef.temp.sav'
;  IF yr EQ 2002 THEN startmo = 8
;  IF yr EQ 2002 THEN endmo = 8
;  IF yr EQ 2005 THEN startmo = 8
;  IF yr EQ 2008 THEN endmo = 10
if yr eq 2010 then restore,mydocs()+'acme07/data/narr/fef.temp.sav'
if yr eq 2010 then startmo = 10

  openw,ffl,mydocs()+'acme07/data/narr/feftemp.progress.txt',/get_lun

  FOR mo = startmo,endmo DO BEGIN
    dim = days_in_month(mo,y=yr)
;    IF yr EQ 2002 THEN dim = 3
;    IF yr EQ 2008 AND mo EQ 10 THEN dim = 28
    FOR day = 1,dim DO BEGIN 
      doystr = string(yr,format='(i4.4)')+string(mo,format='(i2.2)')+string(day,format='(i2.2)')
      doy = dy_to_jd(doystr)
      print,'Year ',yr,' Doy ',doy
      printf,ffl,doystr
      dday = day
      if yr eq 2009 and doy eq 20 then dday++
      if yr eq 2010 and doy eq 8 then dday++
stop
      t = readnarr(yr,mo,dday,vars,lats=lat,lons=lon)
      airt[doy-1,*,*,*] = t.tmp2m.data
      dewpt[doy-1,*,*,*] = t.dpt2m.data
      rh[doy-1,*,*,*] = t.rh2m.data
      soilt[doy-1,*,*,*] = t.tsoil10_40cm.data
      soilq[doy-1,*,*,*] = t.soill0_10cm.data
      u[doy-1,*,*,*] = t.ugrd10m.data
      v[doy-1,*,*,*] = t.vgrd10m.data
      precip[doy-1,*,*,*] = t.apcp.data
      zi[doy-1,*,*,*] = t.hpbl.data
      press[doy-1,*,*,*] = t.presnsfc.data
      sw[doy-1,*,*,*] = t.dswrfsfc.data
      snowd[doy-1,*,*,*] = t.snod.data
      snowc[doy-1,*,*,*] = t.snowc.data
    ENDFOR
    save,airt,dewpt,rh,soilt,soilq,u,v,precip,zi,press,sw,snowd,snowc,filename=mydocs()+'acme07/data/narr/fef.temp.sav'
  ENDFOR

  free_lun,ffl

;read Jan 1 of next year
;  IF yr NE 2008 THEN BEGIN 
    t = readnarr(yr+1,1,1,vars,lats=lat,lons=lon)
    airt[365,*,*,*] = t.tmp2m.data
    dewpt[365,*,*,*] = t.dpt2m.data
    rh[365,*,*,*] = t.rh2m.data
    soilt[365,*,*,*] = t.tsoil10_40cm.data
    soilq[365,*,*,*] = t.soill0_10cm.data
    u[365,*,*,*] = t.ugrd10m.data
    v[365,*,*,*] = t.vgrd10m.data
    precip[365,*,*,*] = t.apcp.data
    zi[365,*,*,*] = t.hpbl.data
    press[365,*,*,*] = t.presnsfc.data
    sw[365,*,*,*] = t.dswrfsfc.data
    snowd[365,*,*,*] = t.snod.data
    snowc[365,*,*,*] = t.snowc.data
;  ENDIF 
    
;transpose vars
  airt = reform(transpose(airt,[3,0,1,2]),2928)
  dewpt = reform(transpose(dewpt,[3,0,1,2]),2928)
  rh = reform(transpose(rh,[3,0,1,2]),2928)
  soilt = reform(transpose(soilt,[3,0,1,2]),2928)
  soilq = reform(transpose(soilq,[3,0,1,2]),2928)
  u = reform(transpose(u,[3,0,1,2]),2928)
  v = reform(transpose(v,[3,0,1,2]),2928)
  precip = reform(transpose(precip,[3,0,1,2]),2928)
  zi = reform(transpose(zi,[3,0,1,2]),2928)
  press = reform(transpose(press,[3,0,1,2]),2928)
  sw = reform(transpose(sw,[3,0,1,2]),2928)
  snowc = reform(transpose(snowc,[3,0,1,2]),2928)
  snowd = reform(transpose(snowd,[3,0,1,2]),2928)

;timezone shift
  airt = airt[2:2921,*,*]-273.15
  dewpt = dewpt[2:2921,*,*]-273.15
  rh = rh[2:2921,*,*]
  soilt = soilt[2:2921,*,*]-273.15
  soilq = soilq[2:2921,*,*]
  u = u[2:2921,*,*]
  v = v[2:2921,*,*]
  precip = precip[2:2921,*,*]
  zi = zi[2:2921,*,*]
  press = press[2:2921,*,*]
  sw = sw[2:2921,*,*]
  snowd = snowd[2:2921,*,*]
  snowc = snowc[2:2921,*,*]

  ws = sqrt(u^2 + v^2)
  vp = 100. * 6.112 * exp ( (17.67 * dewpt) / (243.5 + dewpt) )
  es = 100. * 6.112 * exp ( (17.67 * airt) / (243.5 + airt) )
  ess = 100. * 6.112 * exp ( (17.67 * soilt) / (243.5 + soilt) )
  vpd = (es-vp)>0
  vpds = ess-vp

  tms = 1.0+(findgen(365.*8.)/8.)

  SW_TO_PAR = 2.03 * 10800./1e6  ;* 86400./1e6
  par = sw*SW_TO_PAR

  clim = fltarr(14,365*8.)
  clim2 = fltarr(14,365*8)

;Loc,Year,Day,Hour,Timestep

  clim[0,*] = 0.
  clim[1,*] = yr
  clim[2,*] = (indgen(365*8)/8)+1
  clim[3,*] = (indgen(365*8) MOD 8)*3
  clim[4,*] = 0.125
  clim[5,*] = airt
  clim[6,*] = soilt
  clim[7,*] = sw * SW_TO_PAR 
  clim[8,*] = precip
  clim[9,*] = vpd
  clim[10,*] = vpds
  clim[11,*] = vp
  clim[12,*] = ws
  clim[13,*] = soilq

  clim2[0,*] = 0.
  clim2[1,*] = yr
  clim2[2,*] = (indgen(365*8)/8)+1
  clim2[3,*] = (indgen(365*8) MOD 8)*3
  clim2[4,*] = 0.125
  clim2[5,*] = press/1e2
  clim2[6,*] = rh
  clim2[7,*] = dewpt
  clim2[8,*] = sw
  clim2[9,*] = u
  clim2[10,*] = v
  clim2[11,*] = snowc
  clim2[12,*] = snowd
  clim2[13,*] = zi

  h1 = ['Loc','Year','Day','Hour','Timestep','AirT','SoilT','PAR','Precip','VPD','VPD_Soil','VP','Wind_spd','Soil_moist']
  h2 = ['Loc','Year','Day','Hour','Timestep','Press','RH','DewPt','Shortwave','U','V','Snow_Cov','Snow_D','zi']

  write_xdf,mydocs()+'acme07/data/narr/fef1_'+yrstr+'.xdf',clim,header=h1
  write_xdf,mydocs()+'acme07/data/narr/fef2_'+yrstr+'.xdf',clim2,header=h2

  g = where(~finite(clim),ng)
  IF ng GT 0 THEN clim[g] = -999.0
  g = where(~finite(clim2),ng)
  IF ng GT 0 THEN clim2[g] = -999.0

  form = '(i3," ",i4," ",i3," ",i2," ",f0.3," ",f14.3,f14.3,f14.3,f14.3,f14.3,f14.3,f14.3,f14.3,f14.3)'
  openw,fl,mydocs()+'acme07/data/narr/fef'+yrstr+'.clim',/get_lun
  printf,fl,clim,format=form
  free_lun,fl
  openw,fl,mydocs()+'acme07/data/narr/fef'+yrstr+'.clim2',/get_lun
  printf,fl,clim2,format=form
  free_lun,fl

END