; (3/1/2005)
; ERA40 2D invariant fields (ds118.0) conversion from GRIB to NetCDF
;
; NOTES
; -----
; 1. These fields are on a regular 2.5 deg grid (144x73)
; 2. Field list:
;        Var_name         PN    Name            Long_name
;        --------         --    ----            ---------
;    Z_GDS0_HYBL         129    Z       Surface geopotential
;    LSM_GDS0_SFC        172    LSM     Land-sea mask

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

; --------------- Main block ---------------

begin

; ..Read working directory from environment variable:
;   setenv DS118 ${wdir} (set in ds118.0.csh)
  cla   = getenv("DS118")
  clc   = stringtochar(cla)
  lst   = dimsizes(clc)-2
  if (clc(lst) .eq. "/") then
    lst = lst-1
  end if
  tdir  = chartostring(clc(0:lst))
  system("mkdir -p "+tdir)
  wdir  = tdir+"/"

; ..Read Grib data from MSS to local disk
  locn = "e4inv.fixed"
  system("msrcp -n mss:/DSS/U01487 "+wdir+locn)

; ..Regular grid
  nlat  = 73
  mlon  = 144
  lat   = latGlobeF (nlat, "lat", "latitdue", "degrees_north")
  lon   = lonGlobeF (mlon, "lon", "longitude", "degrees_east")

; ..Open grib file
  grib_in = addfile(wdir+locn+".grb","r")

; ..Create time for HOPS
  time   = (/ 12 /)
  ntim   = dimsizes(time)
  time!0 = "time"
  time@long_name = "initial time"
  time@units     = "hours since 1800-01-01 00:00"

  date   = (/ 19570901 /)
  date!0 = "time"
  date@long_name = "current date as 8 digit integer (YYYYMMDD)"

  datesec   = 3600 * time
  datesec!0 = "time"
  datesec@long_name = "current seconds of current date"
  datesec@units     = "s"

; ..Read two scalar variables: GRIB lats run N->S
  Zg   = grib_in->Z_GDS0_HYBL
  LSMg = grib_in->LSM_GDS0_SFC

; ..Make lat coords S->N
  Zg   = Zg(::-1,:)
  LSMg = LSMg(::-1,:)

; ..Add CAM2-like meta-data
  dimm          = dimsizes(Zg)
  nlat          = dimm(0) 
  mlon          = dimm(1) 
  Z             = new((/1,nlat,mlon/),float)
  LSM           = new((/1,nlat,mlon/),float)
  Z!0           = "time"
  Z!1           = "lat"
  Z!2           = "lon"
  Z&lat         = Zg&g0_lat_2               ; now assign coordinate variables
  Z&lon         = Zg&g0_lon_3               ; to our new, expanded-in-time field
  Z@long_name   = Zg@long_name
  Z@units       = Zg@units
  LSM!0         = "time"
  LSM!1         = "lat"
  LSM!2         = "lon"
  LSM&lat       = LSMg&g0_lat_0
  LSM&lon       = LSMg&g0_lon_1
  LSM@long_name = LSMg@long_name
  LSM@units     = LSMg@units

; ..Copy field data
  Z(0,:,:)   = (/ Zg /)
  LSM(0,:,:) = (/ LSMg /)

; ..Print out summary of variables
  printVarSummary(Z)
  printVarSummary(LSM)

; ..Open output NetCDF file; create parameters and output file
  system("rm -f "+wdir+locn+".nc")
  netcdf_out = addfile(wdir+locn+".nc","c")

; ..Assign file attributes
  nl         = inttochar(10)              ; carriage return
  fAtt       = True
  fAtt@title = "ERA40 GRIB-to-NetCDF"
  fAtt@story = nl+\
    "ERA40 GRIB file: These are invariant 2D fields on a 2.5 degree regular  "+nl+\
    "grid (144x73).  We simply read in the two variables that we want and    "+nl+\
    "write them out to another NetCDF file after changing or removing some   "+nl+\
    "of the meta-data.                                                       "
  fAtt@Conventions   = "None"
  fAtt@source_file   = locn
  fAtt@source_center = "European Center for Medium-Range Weather Forecasts - Reading"
  fAtt@creation_date = systemfunc("date")
  fileattdef( netcdf_out, fAtt )

; ..Predefine coordinate information
  dimNames = (/ "time", "lat", "lon" /)
  dimSizes = (/ ntim  , nlat ,  mlon /)
  dimUnlim = (/ True  , False, False /)
  filedimdef(netcdf_out, dimNames, dimSizes, dimUnlim)

  filevardef   (netcdf_out, "time", typeof(time), "time")
  filevarattdef(netcdf_out, "time", time)

  filevardef   (netcdf_out, "lat" , typeof(lat) , "lat")
  filevarattdef(netcdf_out, "lat" , lat)

  filevardef   (netcdf_out, "lon" , typeof(lon) , "lon")
  filevarattdef(netcdf_out, "lon" , lon)

; ..Predefine variable sizes
  filevardef   (netcdf_out, "date"   , typeof(date)   , getvardims(date))
  filevarattdef(netcdf_out, "date"   , date)
  filevardef   (netcdf_out, "datesec", typeof(datesec), getvardims(datesec))
  filevarattdef(netcdf_out, "datesec", datesec)
  filevardef   (netcdf_out, "Z"      , typeof(Z)      , getvardims(Z))
  filevarattdef(netcdf_out, "Z"      , Z)
  filevardef   (netcdf_out, "LSM"    , typeof(LSM)    , getvardims(LSM))
  filevarattdef(netcdf_out, "LSM"    , LSM)

; ..Write data values to predefined locations
  netcdf_out->time    = (/ time    /)
  netcdf_out->lat     = (/ lat     /)
  netcdf_out->lon     = (/ lon     /)
  netcdf_out->date    = (/ date    /)
  netcdf_out->datesec = (/ datesec /)
  netcdf_out->Z       = (/ Z       /)
  netcdf_out->LSM     = (/ LSM     /)

  print("Normal termination of program ds118.0.ncl")
end