// condense_ccm.i: Condense monthly mean CCM3 data

// Reads monthly mean hybrid level data from CCM3 history files.
// This script reads selected variables from monthly-mean hybrid-level
// CCM3 history files, and interpolates the data to the 17 NCEP
// reanalysis pressure levels:
// 1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10 hPa.
// 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"

// Run name
run_name= "SSTens2";

// Case name (may be different from the run name for name changing runs)
case_name= run_name;

// Starting years for condensed archives
// (the last year on the list need not have data associated with it)
yearlist= [1, 16, 31, 46, 51];

// Directory conventions ("csm"/"ccm3"/"ccm2")
dir_conventions= "ccm3";

// Monthly average file name conventions ("yyyy-mm"/"mm-yy")
mavg_conventions= "mm-yy";

// Input/output user names
//inuser= "MVERTENS";
inuser= "JEFF";
outuser= "SVN";

// Data URL
//data_URL= "http://www.cgd.ucar.edu/csm/exp/b003.html";
data_URL= "http://www.cgd.ucar.edu/gds/jeff/jeff.html";

// Pressure levels for interpolation (hPa)
plevs= double([1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10]);

// 2-D fields to be extracted
var2d= [ "TS","PS","ORO","TREFHT","PRECL","PRECC",
         "SHFLX","LHFLX","QFLX","USTAR",
         "FSNS","FLNS","FSNT","FLNT",
         "FSNSC","FLNSC","FSNTC","FLNTC",
         "TAUX","TAUY","SRFRAD","CLDTOT",
         "BEVAP","BTRAN" ];

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

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

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

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

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

// Input directory name
if (dir_conventions == "csm") {
  indir=  "csm/"+run_name+"/ccm3/hist";
} else {
  indir=  dir_conventions+"/"+run_name+"/hist";
}

// 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
  yearstr= strmid(strnum(10000+startyear,format="%5d"),1,4);
  outfile= "h"+strnum(endyear-startyear+1)+"-"+yearstr+".nc";
  
  outstruc= NULL;
  for (iyear=startyear; iyear <= endyear; iyear++) {
  
    // Year string
    yearstr= strmid(strnum(10000+iyear,format="%5d"),1,4);
  
    // List of filenames for each month
    monfile= array(string,12);
    for (imon=1; imon <= 12; imon++) {
      monstr= strmid(strnum(100+imon,format="%3d"),1,2);
      if (mavg_conventions == "yyyy-mm") {
        monfile(imon)= inprefix + "/" + yearstr+ "-" + monstr;
      } else {
        monfile(imon)= inprefix + "/" + monstr+ "-" + strmid(yearstr,2,2);
      }
    }
  
    // 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
      SLICE= hget(var2d);
  
      // Extract 3D variables
      VARXYZ= hget(var3d);
      VAV= NULL;
      XAV= NULL;
  
      // Compute mass-weighted vertical averages
      for (ivar=1; ivar<=numberof(var_vav); ivar++) {
        islab= strloc(var3d,var_vav(ivar));
        if (islab == 0)
          error, "Cannot compute vertical average of "+var_vav(ivar);

        grow, VAV, hsub(VARXYZ(islab),z="avg",strip="z");
      }
  
      // Interpolate to pressure levels (hPa)
      pvals= hop( hcoord(VARXYZ(1),"z"), "/", 100., units="hPa" );
      VARXYZ= hinterp(VARXYZ, "z", coord=pvals, op="log", grid=plevs);
  
      // 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;

      // Set names for averaged fields
      XAV.name= XAV.name + "_XAV";
      VAV.name= VAV.name + "_VAV";
  
      // Set global attributes
      hset_attr, SLICE, ":data_URL", data_URL;
      hset_attr, XAV,   ":data_URL", data_URL;
      hset_attr, VAV,   ":data_URL", data_URL;
  
      // Save/append to output file
      happend, outstruc, create=outprefix+"/"+outfile, SLICE, XAV, VAV;
  
      // 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