%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  oct_process_coast_CFSR.m
%
% Preprocessing script to remove the land values in CFSR data and 
% extrapolate coastal values inland.
%
% It also removes the redondant values around the Greenwitch Meridian
%
% This is done before using online interpolation
%
%  This file is part of CROCOTOOLS
%
%  CROCOTOOLS is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published
%  by the Free Software Foundation; either version 2 of the License,
%  or (at your option) any later version.
%
%  CROCOTOOLS is distributed in the hope that it will be useful, but
%  WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
%  MA  02111-1307  USA
%
%  Copyright (c) 2012 by Pierrick Penven
%  e-mail:Pierrick.Penven@ird.fr
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear all
close all
%%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%%
%
% Common parameters
%
crocotools_param
%
% Directory of the input CFSR unprocessd data
%
CFSR_dir=[FORC_DATA_DIR,'CFSR_',CROCO_config,'/'];
%
% Directory of the output CFSR processd data
%
CFSR_COAST_dir=[FORC_DATA_DIR,'CFSR_COAST_',CROCO_config,'/'];
%
% Width of the extrapolation band inland (200km is ok for a 1/4deg simulation).
%
distmax=200e3;
%
% Decorrelation scale for the objective analysis [m]
%
ro=2000e3;
%
% CFSR atmospheric variables names
%
vnames={'Land_cover_1land_2sea' ...    % surface land-sea mask [1=land; 0=sea] 
        'Temperature_height_above_ground' ...      % 2 m temp. [k]
        'Downward_Long-Wave_Rad_Flux' ...   % surface downward long wave flux [w/m^2] 
        'Upward_Long-Wave_Rad_Flux_surface' ...   % surface upward long wave flux [w/m^2] 
        'Temperature_surface' ...     % surface temp. [k]   
	'Downward_Short-Wave_Rad_Flux_surface' ...   % surface downward solar radiation flux [w/m^2]
	'Upward_Short-Wave_Rad_Flux_surface' ...   % surface upward solar radiation flux [w/m^2] 
	'Precipitation_rate' ...   % surface precipitation rate [kg/m^2/s] 
	'U-component_of_wind' ...    % 10 m u wind [m/s]
	'V-component_of_wind' ...    % 10 m v wind [m/s]
	'Specific_humidity'};       % 2 m specific humidity [kg/kg]
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of user input  parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
eval(['!mkdir ',CFSR_COAST_dir])
%
% 1 - Process the land sea mask
%
k=1;
disp(['  Processing mask'])
in_file =  [CFSR_dir,char(vnames(k)),'.nc'];
out_file = [CFSR_COAST_dir,char(vnames(k)),'.nc'];
%
ncid = netcdf.open(in_file, 'NC_NOWRITE');
lon=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lon'));
lat=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lat'));
mask=squeeze(netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'Land_cover_1land_2sea')));
netcdf.close(ncid);
%
% 1.1 - Remove the redondant values around greenwitch meridian
% (around lon=0 there is an overlap that could create problems)
%
ibad=find(lon>-0.01 & lon<0.3);
if isempty(ibad)
  ibmin=1;
  ibmax=2;
else
  ibmin=min(ibad)-1;
  ibmax=max(ibad)+1;
end
%
mask=cat(2,mask(:,1:ibmin),mask(:,ibmax:end));
lon=cat(1,lon(1:ibmin),lon(ibmax:end));
[Mm,Lm]=size(mask);
%
% 1.2 - Write the land sea mask down
%
oct_write_NCEP([CFSR_COAST_dir,char(vnames(k)),'.nc'],...
	    char(vnames(k)),lon,lat,0,mask,Yorig)
%
% 2 - Prepare for the extrapolations
%
% 2.1 - Remove the lakes in Africa and South America (not finished)
%
[X,Y]=meshgrid(lon,lat);
mask(X>10 & X<18 & Y>11 & Y<15)=1;
mask(X>16 & X<36 & Y>-15 & Y<10)=1;
mask(X>16 & X<33 & Y>-21 & Y<-12)=1;
mask(X>13 & X<21 & Y>-4   & Y<  15)=1;
mask(X>-72 & X<-69 & Y>9   & Y<  10.5)=1;
mask(X>-55 & X<-52 & Y>-4  & Y<  -2)=1;
mask(X>-69.5 & X<-68 & Y>-17 & Y<-15)=1;
mask(X>-48 & X<-47.5 & Y>-1.5 & Y<-0.7)=1;
%
% 2.2 - Get the number of points to create a "coastal band"
%
[dx,dy]=oct_get_dx(X,Y);  
dxmin=min([min(dx(:)) min(dy(:))]);
npts=ceil(distmax/dxmin);
disp(['Number of points to perform the extrapolation: ',num2str(npts)])
%
% 2.3 - Create the matrices of good points and missing for the extrapolations
%
m2=mask;
for i=1:npts
  m2=oct_hanning(m2);
end
%
mgood=double(m2<1 & m2>0 & mask==0);
mgood=double(m2<1 & m2>0 & mask==0);
mbad=double(m2<1 & m2>0 & mask==1);
mbad=double(m2<1 & m2>0 & mask==1);
disp(['Number of points used for the extrapolation: ',num2str(sum(mgood(:)))])
disp(['Number of points extrapolated: ',num2str(sum(mbad(:)))])
%
% 2.4 - Create the matrix for the extrapolation
%
if 1==1
  disp('Create the OA matrix')
  tic
  coef = oct_oacoef(X(mgood==1),Y(mgood==1),X(mbad==1),Y(mbad==1),ro);
  toc
  save coef.mat coef
else
  disp('Load the OA matrix')
  load coef.mat
end
%
% 3 - Loop on the years
%
for Ye=Ymin:Ymax
  disp(['=========================='])
  disp(['Processing year: ',num2str(Ye)])
  disp(['=========================='])
%
% 3.1 - Loop on the months
%
  if Ye==Ymin
    mo_min=Mmin;
  else
    mo_min=1;
  end
  if Ye==Ymax
    mo_max=Mmax;
  else
    mo_max=12;
  end
  
  for Mo=mo_min:mo_max
    disp(['  Processing month: ',num2str(Mo)])
    disp(['=========================='])
%
% 3.2 - Loop on variable names
%
    for k=2:length(vnames)
      vname=char(vnames(k));
      disp(['=========================='])
      disp(['    VNAME IS ',vname]);
      disp(['=========================='])

      in_file = [CFSR_dir,vname,'_Y',num2str(Ye),'M',num2str(Mo),'.nc'];

      ncid = netcdf.open(in_file, 'NC_NOWRITE');
      time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'time'));
      
      N=length(time);
      var=nan*zeros(N,Mm,Lm);

      for tndx=1:N

        disp(['Processing ',num2str(tndx),' of ',num2str(N)])

        var0=squeeze(netcdf.getVar(ncid, netcdf.inqVarID(ncid, vname)));
        var0=cat(2,var0(:,1:ibmin),var0(:,ibmax:end));
	mvar0=mean(var0(mask==0));
	var0(mask==1)=mvar0;
	var0(mbad==1)=mvar0+coef*(var0(mgood==1)-mvar0);
	
	var(tndx,:,:)=var0;

      end
      
      netcdf.close(ncid);
%
% 3.3 - Write in a file
%
      oct_write_NCEP([CFSR_COAST_dir,vname,'_Y',num2str(Ye),...
	          'M',num2str(Mo),'.nc'],...
	          vname,lon,lat,time,var,Yorig)
%
    end % loop k
  end % end loop month
end % end loop year
%
return