% Generate NCEP Gaussian input fields for cube-sphere,
% bulk-formulae, sea-ice integration: a monthly mean
% climatology based on Charmaine's ncep_rgau.

clear all, close all
nx=192; ny=94; nt=1464;
pn='/host/nemo/data2/odap/ncep_rgau/';

% Surface (10-m) zonal wind velocity in m/s
fn='/RG_u10m'; fout='uwind_192_94_12.bin'; disp(fout)
fld=zeros(nx,ny,12); n=0;
for year=1978:2002,  n=n+1;
  yr=year-100*floor(year/100);
  dc=10*floor(yr/10);
  fin=[pn myint2str(dc) fn myint2str(yr)];
  disp(fin)
  tmp=readbin(fin,[nx ny nt]);
  for mo=1:12
    it=(1:122)+(mo-1)*122;
    fld(:,:,mo)=fld(:,:,mo)+mean(tmp(:,:,it),3);
  end
end
fld=fld/n; writebin(fout,fld);

% Surface (10-m) meridional wind velocity in m/s
fn='/RG_v10m'; fout='vwind_192_94_12.bin'; disp(fout)
fld=zeros(nx,ny,12); n=0;
for year=1978:2002,  n=n+1;
  yr=year-100*floor(year/100);
  dc=10*floor(yr/10);
  fin=[pn myint2str(dc) fn myint2str(yr)];
  disp(fin)
  tmp=readbin(fin,[nx ny nt]);
  for mo=1:12
    it=(1:122)+(mo-1)*122;
    fld(:,:,mo)=fld(:,:,mo)+mean(tmp(:,:,it),3);
  end
end
fld=fld/n; writebin(fout,fld);

% Surface (2-m) air temperature in deg K
fn='/RG_tmp2m'; fout='atemp_192_94_12.bin'; disp(fout)
fld=zeros(nx,ny,12); n=0;
for year=1978:2002,  n=n+1;
  yr=year-100*floor(year/100);
  dc=10*floor(yr/10);
  fin=[pn myint2str(dc) fn myint2str(yr)];
  disp(fin)
  tmp=readbin(fin,[nx ny nt]);
  for mo=1:12
    it=(1:122)+(mo-1)*122;
    fld(:,:,mo)=fld(:,:,mo)+mean(tmp(:,:,it),3);
  end
end
fld=fld/n; writebin(fout,fld);

% Surface (2m) specific humidity in kg/kg
fn='/RG_spfh2m'; fout='aqh_192_94_12.bin'; disp(fout)
fld=zeros(nx,ny,12); n=0;
for year=1978:2002,  n=n+1;
  yr=year-100*floor(year/100);
  dc=10*floor(yr/10);
  fin=[pn myint2str(dc) fn myint2str(yr)];
  disp(fin)
  tmp=readbin(fin,[nx ny nt]);
  for mo=1:12
    it=(1:122)+(mo-1)*122;
    fld(:,:,mo)=fld(:,:,mo)+mean(tmp(:,:,it),3);
  end
end
fld=fld/n; writebin(fout,fld);

% Precipitation in m/s
fn='/RG_rain'; fout='precip_192_94_12.bin'; disp(fout)
fld=zeros(nx,ny,12); n=0;
for year=1978:2002,  n=n+1;
  yr=year-100*floor(year/100);
  dc=10*floor(yr/10);
  fin=[pn myint2str(dc) fn myint2str(yr)];
  disp(fin)
  tmp=readbin(fin,[nx ny nt]);
  for mo=1:12
    it=(1:122)+(mo-1)*122;
    fld(:,:,mo)=fld(:,:,mo)+mean(tmp(:,:,it),3);
  end
end
fld=fld/n; writebin(fout,fld);

% Downward shortwave radiation in W/m^2
fn='/RG_dsw'; fout='swdown_192_94_12.bin'; disp(fout)
fld=zeros(nx,ny,12); n=0;
for year=1978:2002,  n=n+1;
  yr=year-100*floor(year/100);
  dc=10*floor(yr/10);
  fin=[pn myint2str(dc) fn myint2str(yr)];
  disp(fin)
  tmp=readbin(fin,[nx ny nt]);
  for mo=1:12
    it=(1:122)+(mo-1)*122;
    fld(:,:,mo)=fld(:,:,mo)+mean(tmp(:,:,it),3);
  end
end
fld=fld/n; writebin(fout,fld);

% Downward longwave radiation in W/m^2
fn='/RG_dlw'; fout='lwdown_192_94_12.bin'; disp(fout)
fld=zeros(nx,ny,12); n=0;
for year=1978:2002,  n=n+1;
  yr=year-100*floor(year/100);
  dc=10*floor(yr/10);
  fin=[pn myint2str(dc) fn myint2str(yr)];
  disp(fin)
  tmp=readbin(fin,[nx ny nt]);
  for mo=1:12
    it=(1:122)+(mo-1)*122;
    fld(:,:,mo)=fld(:,:,mo)+mean(tmp(:,:,it),3);
  end
end
fld=fld/n; writebin(fout,fld);