%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   THIS PROGRAM : 
%        . Prepare files to Correct the PDF bias during MJJASO
%        . Match GCM PDF at each individual station 
%        . Then use a batch file to creat the script for DOS c:/temp/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Correction prog from Ines correct month by month : mean values : 
% 	1. amount 
%	2. Intensity
%	3. Frequency 

clear;
close all;
stat=['Bakel';'Bambe';'Dakar';'Diour';'Fatic';'Kaola';'Kedou';'kolda';'Lingu';'Louga';
    'Matam';'Mbour';'Nioro';'Podor';'StLou';'Tamba';'Thies';'Velin';'Zigui';'Simen'];
% load /beluga/data/ousmane/data/kergi/station/daily/stat/latlon.txt;
ms=xlsread('/beluga/data/ousmane/data/kergi/station/daily/stat/latlon.xls');
lon=ms(:,1);
lat=ms(:,2);

imon=['Jan';'Feb';'Mar';'Apr';'May';'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec'];
ind0=[1 2 3 4 5 6 9 10 11 12 14 15 17]; % matam + bakel added=11+1
ind1=[7 8 13 16 18 19 20];
iwhole=[1:20];
station=['iwhole';'inorth';'isouth'];
Grille=['Grwhole';'Grnorth';'Grsouth'];
in1=[15 10 9 17 3 12 2 4 5 6]; % missing 6=Klck and 14=Podor
in2=[11 1];
in3=[19 8 13 18];
in4=[7 16 20];
Grwhole=[2 4 5 6 8 9]; % index of grid box covering senegal
Grnorth=[2 4 5];
Grsouth=[6 8 9];
Grid=['All 20 Stations    ';'11 Northern Station';' 9 Southern Station'];
region=['Middle-West 6 ';'North 1       ';'Middle 4      ';'Middle-South 4';'Middle-East 2 ';'South-East 3  '];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        MAY ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load /beluga/data/ousmane/mos_iri/CFS/matlab/new/CFS_St30_30Apr04May_ICs09xyAprOct8108.mat;
cfs0=zeros([4,5,184,28]);
cfs0(1,:,4:184,:)=squeeze(pprate(:,2,2,34:214,:));
cfs0(2,:,4:184,:)=squeeze(pprate(:,2,3,34:214,:));
cfs0(3,:,4:184,:)=squeeze(pprate(:,3,2,34:214,:));
cfs0(4,:,4:184,:)=squeeze(pprate(:,3,3,34:214,:));
% data starts : from 30 Apr=memb1 to 04May=memb5 missing=NaN

load /beluga/data/ousmane/data/kergi/station/daily/data/xamjjaso20stat5008.mat;
xmjjaso=xamjjaso(1:20,32:59,31:214); % for 01May-31Oct 1950-2008
clear xamjjaso;

% Writting GCM files in correct format : new format
for igrid=1:4
    for imemb=1:5
        inn=sprintf('in%d',igrid);
        n=eval(inn);
        for is=1:length(n)
           name=sprintf('%s_G%dM%d.txt',stat(n(is),:),igrid,imemb);
           fid=fopen(name,'wt');
           fprintf(fid,'Longitude Latitude Time Total Precipitation\n');
           fprintf(fid,'degree_east degree_north days since 1948-01-01 00:00:00 m/s\n');
           for iy=1981:2008
               xval=ones(1,366)*-99.;
        if((rem(iy,4)==0 & rem(iy,100)~=0) |(rem(iy,400)==0)) % leap years
                   xval(1,122:305)=cfs0(igrid,imemb,1:184,iy-1980)/(86400*1000); % from mm/day -> m/s
                   for id=1:366
                nn=sprintf('datestr(datenum(''1.1.%d'',''dd.mm.yyyy'')+id-2,24)',iy);
                aa=eval(nn);
                dad=sscanf(aa,'%2d/%2d/%4d');
                nn=sprintf('datestr(datenum(''1.1.%d'',''dd.mm.yyyy'')+id-1,24)',iy);
                aa=eval(nn);
                dad1=sscanf(aa,'%2d/%2d/%4d');
                fprintf(fid,'%7.4fW %7.4fN 2100 %d %s %d - 2100 %d %s %d %E\n',-lon(n(is)),lat(n(is)),dad(1),num2mon(dad(2)),dad(3),dad1(1),num2mon(dad1(2)),dad1(3),xval(1,id));
                   end
	else
	        xval(1,121:304)=cfs0(igrid,imemb,1:184,iy-1980)/(86400*1000); % from mm/day -> m/s
                   for id=1:365
                nn=sprintf('datestr(datenum(''1.1.%d'',''dd.mm.yyyy'')+id-2,24)',iy);
                aa=eval(nn);
                dad=sscanf(aa,'%2d/%2d/%4d');
                nn=sprintf('datestr(datenum(''1.1.%d'',''dd.mm.yyyy'')+id-1,24)',iy);
                aa=eval(nn);
                dad1=sscanf(aa,'%2d/%2d/%4d');
                fprintf(fid,'%7.4fW %7.4fN 2100 %d %s %d - 2100 %d %s %d %E\n',-lon(n(is)),lat(n(is)),dad(1),num2mon(dad(2)),dad(3),dad1(1),num2mon(dad1(2)),iy,xval(1,id));
                   end
	end
           end
           fclose(fid);
        end
    end
end

% Now station data in format
for is=1:20
    name=sprintf('o_%s.wtd',stat(is,:));
    fid=fopen(name,'wt');
    fprintf(fid,'@  DATE  SRAD  TMAX  TMIN  RAIN\n');    
    for iy=1981:2008
               xval=ones(1,366)*-99.;
        if((rem(iy,4)==0 & rem(iy,100)~=0) |(rem(iy,400)==0)) % leap years
                   xval(1,122:305)=xmjjaso(is,iy-1980,1:184); % in mm/day
                   for id=1:366                
                nn=sprintf('datestr(datenum(''1.1.%d'',''dd.mm.yyyy'')+id-1,24)',iy);
                aa=eval(nn);
                dad=sscanf(aa,'%2d/%2d/%4d');
                fprintf(fid,'%d%03d%6.1f%6.1f%6.1f%6.1f\n',iy,id,-99.,-99.,-99.,xval(1,id));
                   end
        else
                xval(1,121:304)=xmjjaso(is,iy-1980,1:184); % in mm/day
                   for id=1:365
                nn=sprintf('datestr(datenum(''1.1.%d'',''dd.mm.yyyy'')+id-1,24)',iy);
                aa=eval(nn);
                dad=sscanf(aa,'%2d/%2d/%4d');
                fprintf(fid,'%d%03d%6.1f%6.1f%6.1f%6.1f\n',iy,id,-99.,-99.,-99.,xval(1,id));
                   end
        end
    end
    fclose(fid);
end

% preparing a bat file for DOS program correction
fid=fopen('biasPDF.bat','wt');
fprintf(fid,'@Echo off \n');
fprintf(fid,'Echo Performing GCM Bias correction \n');
for igrid=1:4
    for imemb=1:5
        inn=sprintf('in%d',igrid);
        n=eval(inn);
        for is=1:length(n)
           fprintf(fid,'GCM_Bias_Corr2 o_%s.wtd %s_G%dM%d.txt 1981 2009 121 304 0 0\n',stat(n(is),:),stat(n(is),:),igrid,imemb);
           fprintf(fid,'RENAME %s_cd.pcp c_%s_G%dM%d.txt\n',stat(n(is),:),stat(n(is),:),igrid,imemb);
        end
    end
end
fclose(fid);

% prepare to do the correction under linux
% for igrid=1:1
%     for imemb=1:1
%         inn=sprintf('in%d',igrid);
%         n=eval(inn);
%         for is=1:size(n,2)
%            name=sprintf('%s_G%dM%d.txt',stat(n(is),:),imemb,igrid);:q

%            station=sprintf('%s.wtd',stat(is,:));
%            yeap=sprintf('dos ''GCM_Bias_Corr2 %s.wtd %s 1981 2007 121 304 0''',stat(n(is),:),name)
%            eval(yeap);
%         end
%     end
% end
