// reanalhist.i: Condense monthly mean NCEP reanalysis data

// Reads monthly mean pressure level data from the NCEP Reanalysis
// (http://www.cgd.ucar.edu/cas/catalog/nmc/rean/press/means.html).
// The data are on a T63 gaussian grid (original data were on a 2.5deg grid),
// in CCM history tape format on 17 pressure levels:
// 1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10 hPa.
// This script reads selected variables and truncates the spectral resolution
// to T42. Also, data are sliced at selected pressure levels for different
// variables, and zonal averages are also computed. The condensed data set on
// the T42 gaussian grid (128x64) is saved as a hyperslab netCDF file.

#include "~/cgdyorick/shtran.i"

// Case name
case_name= "REANAL";

// Starting years for condensed archives
// (the last year on the list need not have data associated with it)
yearlist= [87, 89];

// Case title
case_title= "NCEP reanalysis";

// Data URL
data_URL= "http://www.cgd.ucar.edu/cas/catalog/nmc/rean/press/means.html";

// 2-D fields to be extracted
var2d= ["PS"];

// 3-D fields to be extracted
var3d= ["ZA", "Q", "T"];

// Slices of 3-D fields at selected levels (hPa)
var_slice= ["ZA", "Q", "T", "T"];
var_lev=   [ 500, 850, 850, 300];

// Zonal averages of selected 3-D fields
var_xav= ["Q", "T"];

// Vector 3-D fields to be extracted
vec3d= [["U","V"]];

// Slices of vector 3-D fields at selected levels (hPa)
vec_slice= ["U", "U", "U", "V", "V", "V"];
vec_lev=   [850, 300, 200, 850, 300, 200];

// Zonal averages of selected vector 3-D fields
vec_xav= ["U", "V"];

// Retention period in days (for MSS output files)
rtpd= 1100;

// Input/output user names
inuser= "CA";
outuser= "SVN";

// Input directory name
indir= "NMC/PRESS";

// Output directory name
outdir= "csm/"+case_name+"/ccm3/hist";

// Temporary directory
tmpdir= get_env("TMPDIR");
if (strlen(tmpdir) == 0) tmpdir=".";
write, "TMPDIR="+tmpdir;

// Directory prefixes
inprefix=  inuser+"/"+indir;
outprefix= outuser+"/"+outdir;

// Temporary disk directories for input/output files
oscommand, "mkdir -p "+tmpdir+"/"+inprefix;
oscommand, "mkdir -p "+tmpdir+"/"+outprefix;

for (iperiod=1; iperiod <= numberof(yearlist)-1; iperiod++) {
  
  // Staring/ending years
  startyear= yearlist(iperiod);
  endyear= yearlist(iperiod+1)-1;
  
  // Output file name
  outfile= "h"+strnum(endyear-startyear+1)+"-19"+strnum(startyear)+".nc";
  
  outstruc= NULL;
  for (iyear=startyear; iyear <= endyear; iyear++) {
  
    // Year string
    yearstr= strmid(strnum(10000+iyear,format="%5d"),3,2);
  
    // List of filenames for each month
    monfile= array(string,12);
    for (imon=1; imon <= 12; imon++) {
      monfile(imon)= inprefix+"/"+"T63"+yearstr+strmid(strnum(100+imon,format="%3d"),1,2)+"M"; }
  
    // cd to temporary working directory
    cd, tmpdir;
  
    // Pre-read history files from MSS
    oscommand, "msreadlist -nc "+strcombine(monfile(*)," ");
  
    for (imon=1; imon <= 12; imon++) {
      // Open monthly file after converting from CCM history to netCDF format
      // (Delete CCM history file after conversion, and netCDF file at end)
  
      hopen, monfile(imon)+".nc", case_name=case_name,newtime="date",scratch=1;
  
      // Extract 2-D scalar variables on T63 grid
      SLICE= hget(var2d);
  
      // Truncate to T42 grid
      SLICE= hshop(SLICE, m1=42, n1=42, precision="float");

      // Extract 3D variables on T63 grid
      VARXYZ= hget(var3d);
      XAV= NULL;

      // Truncate to T42 grid
      VARXYZ= hshop(VARXYZ, m1=42, n1=42, precision="float");

      iqslab= strloc(var3d,"Q");
      if (iqslab > 0) {
        // Q field; re-introduce missing values
        VARXYZ(iqslab)= hmask(VARXYZ(iqslab), hop(VARXYZ(iqslab),"<",1.0e33));
      }

      // Extract slices
      for (ivar=1; ivar<=numberof(var_slice); ivar++) {
        islab= strloc(var3d,var_slice(ivar));
        if (islab == 0)
           error,"Cannot extract slice "+var_slice(ivar)+strnum(var_lev(ivar));

        tem_slab= hsub(VARXYZ(islab),z=var_lev(ivar));
        tem_slab.name= tem_slab.name + strnum(var_lev(ivar));
        grow, SLICE, tem_slab;
      }

      // Compute zonal averages
      for (ivar=1; ivar<=numberof(var_xav); ivar++) {
        islab= strloc(var3d,var_xav(ivar));
        if (islab == 0)
          error, "Cannot compute X-average of "+var_xav(ivar);

        grow, XAV, hsub(VARXYZ(islab),x="avg");
      }

      // Delete 3-D variables
      VARXYZ= NULL;

      // Extract vector 3D variables on T63 grid
      VECXYZ= hget(vec3d);

      // Truncate to T42 grid
      VECXYZ= hshop(VECXYZ, vec=1, m1=42, n1=42, precision="float");

      // Extract slices
      for (ivar=1; ivar<=numberof(vec_slice); ivar++) {
        islab= strloc(vec3d(*),vec_slice(ivar));
        if (islab == 0)
           error,"Cannot extract slice "+vec_slice(ivar)+strnum(vec_lev(ivar));

        tem_slab= hsub(VECXYZ(islab),z=vec_lev(ivar));
        tem_slab.name= tem_slab.name + strnum(vec_lev(ivar));
        grow, SLICE, tem_slab;
      }

      // Compute zonal averages
      for (ivar=1; ivar<=numberof(vec_xav); ivar++) {
        islab= strloc(vec3d(*),vec_xav(ivar));
        if (islab == 0)
          error, "Cannot compute X-average of "+vec_xav(ivar);

        grow, XAV, hsub(VECXYZ(islab),x="avg");
      }

      // Delete 3-D vectors
      VECXYZ= NULL;

      // Set names for averaged fields
      XAV.name= XAV.name + "_XAV";
  
      // Set global attributes
      hset_attr, SLICE, ":case_title", case_title;
      hset_attr, XAV,   ":case_title", case_title;
      hset_attr, SLICE, ":data_URL",   data_URL;
      hset_attr, XAV,   ":data_URL",   data_URL;
  
      // Save/append to output file
      happend, outstruc, create=outprefix+"/"+outfile, SLICE, XAV;
  
      // Close input file
      hclose;
    }
  
  }

  // Close output file and save to mass store
  hclose, outstruc, command="cd "+tmpdir+"/"+outuser+"; mssput -rtpd "+strnum(rtpd)+" "+outdir+"/"+outfile;
  
  write, "Output file /"+outprefix+"/"+outfile+" saved to mass store";

}
  
quit