clear
close all;tic;


fname_rain_cmap='/tigress/GEOCLIM/LRGROUP/shared_data/precipitation/CMAP_2020/precip.mon.mean_enhanced_monthly.nc';
fname_rain_gpcp='/tigress/GEOCLIM/LRGROUP/shared_data/precipitation/GPCPv2.3/precip.mon.mean.nc';
fname_sst='/tigress/GEOCLIM/LRGROUP/shared_data/sst_OISSTmonthly/allsst.avhrr-only-v2.1_1981_2020_monthly.nc';
fname_sst_ano='/tigress/GEOCLIM/LRGROUP/shared_data/sst_OISSTmonthly/allsst.avhrr-only-v2.1_1981_2020_monthly_detrend_deseason.nc';

fn_out=['obs_sst_rain.nc'];
fco2_region=[0 360 -90 90];
trange=[datenum(1982,1,1) datenum(2019,12,31,23,59,59)];

toc;disp(['process file: ', fname_rain_cmap]); tic
[rain,lon_rain,lat_rain,time_rain]=read_ncall('fn',fname_rain_cmap,'vn','precip','xn','lon','yn','lat','tn','time','tb',datenum(1800,1,1),'ts',1/24,'xr',fco2_region(1:2),'yr',fco2_region(3:4),'ry','t','tr',trange);
rain_ano=compute_anomaly_bandpass(rain,time_rain,12*15,13);
% interpolate to a uniform grid
[lon_outgrid,lat_outgrid]=meshgrid([0.5:1:359.5],[-89.5:1:89.5]);
lm=length(time_rain); [jm,im]=size(lon_outgrid);
rain_outgrid=nan(jm,im,lm);         
rain_ano_outgrid=nan(jm,im,lm);
for l=1:lm
    rain_outgrid(:,:,l)=interp2(lon_rain,lat_rain,squeeze(rain(:,:,l)),lon_outgrid,lat_outgrid);
    rain_ano_outgrid(:,:,l)=interp2(lon_rain,lat_rain,squeeze(rain_ano(:,:,l)),lon_outgrid,lat_outgrid);
end
%write to nc file
[var_all]=load_var_all(1);
create_mom6_1var3d(fn_out,'lon','lat','time','time_bnds',var_all(1,:),lon_outgrid(1,:),lat_outgrid(:,1),time_rain-datenum(1900,1,1),[],permute(rain_outgrid,[2,1,3]));
create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(2,:),lon_outgrid(1,:),lat_outgrid(:,1),time_rain-datenum(1900,1,1),[],permute(rain_ano_outgrid,[2,1,3]));

toc;disp(['process file: ', fname_rain_gpcp]); tic
%process the second file
[rain,lon_rain,lat_rain,time_rain]=read_ncall('fn',fname_rain_gpcp,'vn','precip','xn','lon','yn','lat','tn','time','tb',datenum(1800,1,1),'ts',1,'xr',fco2_region(1:2),'yr',fco2_region(3:4),'tr',trange);
rain_ano=compute_anomaly_bandpass(rain,time_rain,12*15,13);
% interpolate to a uniform grid
lm=length(time_rain); [jm,im]=size(lon_outgrid);
rain_outgrid=nan(jm,im,lm);
rain_ano_outgrid=nan(jm,im,lm);
for l=1:lm
    rain_outgrid(:,:,l)=interp2(lon_rain,lat_rain,squeeze(rain(:,:,l)),lon_outgrid,lat_outgrid);
    rain_ano_outgrid(:,:,l)=interp2(lon_rain,lat_rain,squeeze(rain_ano(:,:,l)),lon_outgrid,lat_outgrid);
end
%write to nc file
create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(3,:),lon_outgrid(1,:),lat_outgrid(:,1),time_rain-datenum(1900,1,1),[],permute(rain_outgrid,[2,1,3]));
create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(4,:),lon_outgrid(1,:),lat_outgrid(:,1),time_rain-datenum(1900,1,1),[],permute(rain_ano_outgrid,[2,1,3]));

%process the third file
toc;disp(['process file new: ', fname_sst]); tic
[sst,lon_sst,lat_sst,time_sst]=read_ncall('fn',fname_sst,'vn','sst','xn','lon','yn','lat','tn','time','tb',datenum(1978,1,1,12,0,0),'ts',1,'xr',fco2_region(1:2),'yr',fco2_region(3:4),'tr',trange);
sst=squeeze(sst);
%sst_ano=compute_anomaly_bandpass(sst,time_sst,12*15,13);
[sst_ano,lon_sst,lat_sst,time_sst]=read_ncall('fn',fname_sst_ano,'vn','sst','xn','lon','yn','lat','tn','time','tb',datenum(1978,1,1,12,0,0),'ts',1,'xr',fco2_region(1:2),'yr',fco2_region(3:4),'tr',trange);
% interpolate to a uniform grid
lm=length(time_sst); [jm,im]=size(lon_outgrid);
sst_outgrid=nan(jm,im,lm);
sst_ano_outgrid=nan(jm,im,lm);
for l=1:lm
    sst_outgrid(:,:,l)=interp2(lon_sst,lat_sst,squeeze(sst(:,:,l)),lon_outgrid,lat_outgrid);
    sst_ano_outgrid(:,:,l)=interp2(lon_sst,lat_sst,squeeze(sst_ano(:,:,l)),lon_outgrid,lat_outgrid);
end
%write to nc file
create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(5,:),lon_outgrid(1,:),lat_outgrid(:,1),time_sst-datenum(1900,1,1),[],permute(sst_outgrid,[2,1,3]));
create_mom6_1var3d_append(fn_out,'lon','lat','time','time_bnds',var_all(6,:),lon_outgrid(1,:),lat_outgrid(:,1),time_sst-datenum(1900,1,1),[],permute(sst_ano_outgrid,[2,1,3]));

ncwriteatt(fn_out,'/','data description',['this is regridded data from CMAP and GPCP rain and OISST ']);


function [var_all]=load_var_all(load_ind)

   if load_ind==1;
      %write fco2, fco2_ano, dpco2, dpco2_ano, pco2, pco2_ano, kw, kw_ano, fco2_kw, fco2_dpco2, fco2_ct, fco2_resi
      var_all{1,1}='rain_cmap';          var_all{1,2}='mm/d';         var_all{1,3}=['Average Monthly Rate of Precipitation from CMAP'];
      var_all{2,1}='rain_ano_cmap';      var_all{2,2}='mm/d';         var_all{2,3}=['Monthly Anomaly of Precipitation from CMAP'];
      var_all{3,1}='rain_gpcp';          var_all{3,2}='mm/d';         var_all{3,3}=['Average Monthly Rate of Precipitation from CMAP'];
      var_all{4,1}='rain_ano_gpcp';      var_all{4,2}='mm/d';         var_all{4,3}=['Monthly Anomaly of Precipitation from CMAP'];
      var_all{5,1}='sst';                var_all{5,2}='degree';       var_all{5,3}=['Monthly sea surface temperature from OISST avhrr only'];
      var_all{6,1}='sst_ano';            var_all{6,2}='degree';       var_all{6,3}=['Monthly Anomaly of SST from OISST avhrr only'];
   end
end





