; (3/1/2005)
; Convert ERA40 3D pressure level variables in GRIB (ds120.1) to
; multiple 2D variables in 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_ISBL_123     129      Z     Geopotential
;    T_GDS0_ISBL_123     130      T     Temperature
;    U_GDS0_ISBL_123     131      U     Zonal wind component
;    V_GDS0_ISBL_123     132      V     Meridional wind component
;    Q_GDS0_ISBL_123     133      Q     Specific humidity

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

; --------------- Subroutines ---------------

undef ("clevi")
function clevi(lev:integer, vlev:integer)

; Return array of indices where each vlev occurs in lev

begin
  mlev = dimsizes(lev)
  nlev = dimsizes(vlev)
  iclev = new((/nlev/),integer)
  do ilev = 0,nlev-1
    j = 0
    do while (vlev(ilev) .ne. lev(j))
      j = j+1
      if (j .gt. mlev) then
        print("ERROR EXIT -- Requested level: "+vlev(ilev)+" not found")
        print("Available levels:")
        print(lev)
        exit
      end if
    end do
    iclev(ilev) = j
  end do
  return iclev
end

undef ("f4Dto3D")
function f4Dto3D(f:file, varName:string, klvl:integer,\
                 ntStrt:integer, ntLast:integer, Stride:integer)

; Return a single specified level (klvl) from a multiple level variable

begin

  dimSizes = filevardimsizes(f,varName)     ; get dimension sizes
  ntim     = (ntLast-ntStrt+Stride)/Stride  ; # output time slices
  nlat     = dimSizes(2)                    ; # output latitudes
  nlon     = dimSizes(3)                    ; # output longitudes

  temp  = new((/nlat,nlon/), float)
  temp  = 0.0
  
  var   = new((/ntim,nlat,nlon/), float)
  var!0 = "time"
  var!1 = "lat"
  var!2 = "lon"

  NT = -1
  do nt=ntStrt,ntLast,Stride
    NT = NT+1 
    temp = f->$varName$(nt,klvl,:,:)
    var(NT,:,:) = (/temp/)                  ; original grid
  end do

  var@long_name = temp@long_name
  var@units     = temp@units

  return(var)
end

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

begin

; ..Specify needed levels for each variable
  Vlev = (/1000,925,850,300,200,100,50,10/)
  Qlev = (/1000,850,700,500,300/)
  Tlev = (/1000,850,700,500,300,200/)
  Zlev = (/1000,850,500,300,200,100,50/)

; ..Specify index of starting and ending times and stride
;   (0 --> 11 with Stride 1 processes whole year of monthly means)
  ntStrt = 0
  ntLast = 11
  Stride = 1

; ..Read starting year, number of years and working directory from
;   environment variable:  setenv DS120 "${istrt};${nyrs};${wdir}"
  cla = getenv("DS120")
  clc = stringtochar(cla)
  lst = dimsizes(clc)-2
  if (clc(lst) .eq. "/") then
    lst = lst-1
  end if
  i = 0
  do while (clc(i) .ne. ";")
    i = i+1
  end do
  ie = i-1
  i = i+1
  do while (clc(i) .ne. ";")
    i = i+1
  end do
  jb = ie+2
  je = i-1
  ke = je+2
  istrt = stringtointeger(chartostring(clc(0:ie)))
  nfils = stringtointeger(chartostring(clc(jb:je)))
  iend  = istrt+nfils-1
  stp = False
  if (istrt .lt. 1958 .or. istrt .gt. 2001) then
    print("START YEAR: "+istrt+" Must have 1958 <= start year <= 2001")
    stp = True
  end if
  if (iend .lt. 1958 .or. iend .gt. 2001) then
    print("END YEAR  : "+iend +" Must have 1958 <=  end year  <= 2001")
    stp = True
  end if
  if (stp) then
    exit
  end if
  tdir  = chartostring(clc(ke:lst))
  system ("mkdir -p "+tdir)
  wdir  = tdir+"/"

; ..Create lists of MSS and local file names
  mssnV = new((/nfils/),string)
  mssnS = new((/nfils/),string)
  locnV = new((/nfils/),string)
  locnS = new((/nfils/),string)
  locno = new((/nfils/),string)
  idate = new((/nfils/),integer)
  i = -1
  do iyr = istrt,istrt+nfils-1
    i = i+1
    mssnV(i)  = "/DSS/U60"+(iyr-1517)
    mssnS(i)  = "/DSS/U65"+(iyr-1565)
    locnV(i)  = "e4moda.dssuv.pl."+iyr
    locnS(i)  = "e4moda.pl.nouv."+iyr
    locno(i)  = "e4moda.pl."+iyr
    idate(i) = iyr*10000+101
  end do

; ..Read all files to local disk before proceeding with processing
  i = -1
  do iyr = istrt,istrt+nfils-1
    i = i+1
    system ("msrcp -n mss:"+mssnV(i)+" "+wdir+locnV(i))
    system ("msrcp -n mss:"+mssnS(i)+" "+wdir+locnS(i))
  end do

; ..Set the regular grid
  nlat  = 73
  nlon  = 144
  lat   = latGlobeF (nlat, "lat", "latitdue", "degrees_north")
  lon   = lonGlobeF (nlon, "lon", "longitude", "degrees_east")

; ..Loop over all files
  do nf=0,nfils-1
    print("==========> nf="+nf+" <==========")
    wcStrt = systemfunc("date")             ; wall clock start

; ..Open grib files
    gribV = addfile(wdir+locnV(nf)+".grb","r")
    gribS = addfile(wdir+locnS(nf)+".grb","r")
    wallClockElapseTime(wcStrt, "Read/open files: "+locnV(nf)+" "+locnS(nf), 0)

; ..Read array of available pressure levels
    lev   = gribV->lv_ISBL1
    levS  = gribS->lv_ISBL1
    do ilev = 0, dimsizes(lev)-1
      if (lev(ilev) .ne. levS(ilev)) then
        print("Levels on two input files do not match")
        exit
      end if
    end do

; ..Calculate indices for needed levels
    Vilev = clevi(lev, Vlev)
    Qilev = clevi(lev, Qlev)
    Tilev = clevi(lev, Tlev)
    Zilev = clevi(lev, Zlev)

; ..Read and create times
    stime = gribV->initial_time0(ntStrt:ntLast)
    time  = gribV->initial_time0_hours(ntStrt:ntLast)
    ntim  = dimsizes(time)

    DATE  = doubletointeger(gribV->initial_time0_encoded(ntStrt:ntLast))
    date  = DATE/100
    date@long_name = "current date as 8 digit integer (YYYYMMDD)"
    date!0= "time"
    print(date)
    if (date(ntStrt) .eq. idate(nf)) then
      print("Date in input file matches specified date")
    else
      print("STOP -- date in input file does not match specified date")
      exit()
    end if

    hh    = (DATE-date*100)
    datesec   = hh*3600
    datesec!0 = "time"
    datesec@long_name = "current seconds of current date"
    datesec@units     = "s"

; ..Print some things out for first file only
    if (nf.eq.0) then
      print(Vilev)
      print(Qilev)
      print(Tilev)
      print(Zilev)
      print(stime)
      print(time)
      print(DATE)
      print(date)
      print(datesec)
    end if

; ..Remove previous output file if it exists and create a new one
    system ("rm -f "+wdir+locno(nf)+".nc")
    netcdf_out = addfile(wdir+locno(nf)+".nc","c")

; ..Assign file attributes
    wcDef      = systemfunc("date")
    nl         = inttochar(10)              ; carriage return
    fAtt       = True
    fAtt@title = "ERA40 GRIB-to-NetCDF"
    fAtt@story = nl+\
      "ERA40 GRIB file: These are monthly mean 2D surface variables on a 2.5   "+nl+\
      "degree regular grid (144x73).  We simply read in the 19 variables that  "+nl+\
      "we want and write them out to another NetCDF file after changing or     "+nl+\
      "removing some of the meta-data.                                         "+nl
    fAtt@Conventions   = "None"
    fAtt@source_file   = locnV(nf)+" "+locnS(nf)
    fAtt@source_center = "European Center for Medium-Range Weather Forecasts - Reading"
    fAtt@creation_date = systemfunc("date")
    fileattdef( netcdf_out, fAtt )
    dimNames = (/ "time", "lat", "lon" /)
    dimSizes = (/ ntim  , nlat ,  nlon /)
    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)
    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) 

; ..Read in an example of each variable to set the size of the output
;   NetCDF file
    varU = f4Dto3D(gribV, "U_GDS0_ISBL_123", Vilev(0), ntStrt, ntLast, Stride)
    varV = f4Dto3D(gribV, "V_GDS0_ISBL_123", Vilev(0), ntStrt, ntLast, Stride)
    do ilev = 0, dimsizes(Vlev)-1
      filevardef   (netcdf_out, "U"+Vlev(ilev), typeof(varU), getvardims(varU))
      filevarattdef(netcdf_out, "U"+Vlev(ilev), varU)
      filevardef   (netcdf_out, "V"+Vlev(ilev), typeof(varV), getvardims(varV))
      filevarattdef(netcdf_out, "V"+Vlev(ilev), varV)    
    end do 
    delete(varU)
    delete(varV)

    varZ = f4Dto3D(gribS, "Z_GDS0_ISBL_123", Zilev(0), ntStrt, ntLast, Stride)
    do ilev = 0, dimsizes(Zlev)-1
      filevardef   (netcdf_out, "Z"+Zlev(ilev), typeof(varZ), getvardims(varZ))
      filevarattdef(netcdf_out, "Z"+Zlev(ilev), varZ) 
    end do 
    delete(varZ)

    varT = f4Dto3D(gribS, "T_GDS0_ISBL_123", Tilev(0), ntStrt, ntLast, Stride)
    do ilev = 0, dimsizes(Tlev)-1
      filevardef   (netcdf_out, "T"+Tlev(ilev), typeof(varT), getvardims(varT))
      filevarattdef(netcdf_out, "T"+Tlev(ilev), varT) 
    end do 
    delete(varT)

    varQ = f4Dto3D(gribS, "Q_GDS0_ISBL_123", Qilev(0), ntStrt, ntLast, Stride)
    do ilev = 0, dimsizes(Qlev)-1
      filevardef   (netcdf_out, "Q"+Qlev(ilev), typeof(varQ), getvardims(varQ))
      filevarattdef(netcdf_out, "Q"+Qlev(ilev), varQ) 
    end do 
    delete(varQ)

    netcdf_out->time    = (/ time    /)
    netcdf_out->lat     = (/ lat     /)
    netcdf_out->lon     = (/ lon     /)
    netcdf_out->date    = (/ date    /)
    netcdf_out->datesec = (/ datesec /)
    wallClockElapseTime(wcDef, "Definition phase took ", 0)

; ..Read standard scalar variables: GRIB lats run N->S
    wcWrite = systemfunc("date")
    do ilev = 0, dimsizes(Vilev)-1
      varU  = f4Dto3D(gribV, "U_GDS0_ISBL_123", Vilev(ilev), ntStrt, ntLast, Stride)
      varV  = f4Dto3D(gribV, "V_GDS0_ISBL_123", Vilev(ilev), ntStrt, ntLast, Stride)
      varU  = varU(:,::-1,:)
      varV  = varV(:,::-1,:)
      VAR1  = "U"+Vlev(ilev)
      VAR2  = "V"+Vlev(ilev)
      netcdf_out->$VAR1$ = (/ varU /)
      netcdf_out->$VAR2$ = (/ varV /)  
      delete(varU)
      delete(varV)   
    end do

    do ilev = 0, dimsizes(Zilev)-1
       varZ = f4Dto3D(gribS, "Z_GDS0_ISBL_123", Zilev(ilev), ntStrt, ntLast, Stride)
       varZ = varZ(:,::-1,:)
       VAR  = "Z"+Zlev(ilev)
       netcdf_out->$VAR$ = (/ varZ /)
       delete(varZ)
    end do
    do ilev = 0, dimsizes(Qilev)-1
       varQ = f4Dto3D(gribS, "Q_GDS0_ISBL_123", Qilev(ilev), ntStrt, ntLast, Stride)
       varQ = varQ(:,::-1,:)
       VAR  = "Q"+Qlev(ilev)
       netcdf_out->$VAR$ = (/ varQ /)
       delete(varQ)
    end do
    do ilev = 0, dimsizes(Tilev)-1
       varT = f4Dto3D(gribS, "T_GDS0_ISBL_123", Tilev(ilev), ntStrt, ntLast, Stride)
       varT = varT(:,::-1,:)
       VAR  = "T"+Tlev(ilev)
       netcdf_out->$VAR$ = (/ varT /)
       delete(varT)
    end do                    
    wallClockElapseTime(wcWrite, "Write phase took ", 0)

; ..Clean up variables
    delete(stime)
    delete(time)
    delete(date)
    delete(datesec)
  end do                                    ; loop on files

  print("Normal termination of program ds120.1.ncl")
end