%Name: getGPMv03B_Ku_holes_m2_cmed_flagEcho
%Purpose: Process raw GPM data (2Ku)
%Latest update: 
%By: Ulrike Romatschke (Original version by Darren Wilton)
%Updated and Modified by MD.Zuluaga (2011/09/26)
%Updated and Modified by S.Brodzik for GPM (2015/03/02)

%Input: 2Ku GPM products
%The input files are the pre-proceesed version of original files by the
%extract_2Ku_getGPM_all.pro IDL code that is needed to be run first!!!!

% Modified 12/09/11 to include reading and interpolation of the
% near surface rain variable as was done with the rain Type
% variable

% Modified 01/10/12 to use only 2A25 modified files to extract data beacuse 
% now the version 7 includes the 2A23-rainType infromation. You
% have to make sure the 2A25 file does contain this info

% Modified v72 10/10/12 to account for extracting the shallow
% convection type of rain from the rainType variable

% Modified v73 05/06/13 Including the SLH latent heating product
%              04/13/13 Read Files with coordinates to avoid in
%	       de-compressing HDF files without data in the selected region
	       
% Modified v74 03/31/15 Output 'shallowRainTypes', 'widthBB', 'heightBB'
	       
%Output: 18 ascii files
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridColumns.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridHeightBB.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridLatitude.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridLongitude.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridNearSurfPhase.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridNearSurfPhaseRaw.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridNearSurfRain.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridRainType.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridRainTypeRaw.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridRows.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridShallowRainType.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.gridWidthBB.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.hdfFileName.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.latitude.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.longitude.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.noSwath.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.numSwaths.txt
%GPM.<orbit>.<YYYYMMDD>.<hhmmss>.reflectivity.txt

%Other codes/functions used: dateSorterGPM.m datestamper.m datestring.m
%   getHDFdate.m get2Dinterp.m getNearSurfRain.m getRaintype.m 
%   loadHDF_new.m loadHDFvdata.m 

%This code runs for 176 levels and takes 6 digit orbits as input. It
%takes whole orbits as input as well as pre cut orbits.

%This code is the first step to process the raw GPM/2Ku data. It cuts out the
%lat long boxes and corrects the geolocation.

%This code does not work for processing whole orbits!!! You must specify
%the region of interest and this cant be the whole globe!

addpath /home/disk/bob/gpm/code/
clear all;
close all;

t=clock;
datestamper(t);  % this Prints the Actual Time
clear t ans

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%!!!!!!!!!!!!!!!!!!If your area extends over the 180 meridian set over180=1. Default is 0. The
%negative value is MinLongEast, the positive is MaxLongEast. The resulting
%output longitude will be positive for all values, e.g. -175 becomes +185.!!!!!!!!!!!!!
over180=0; %default is zero!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
minus40=0; %default is zero!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

%Number of levels in GPM-2Ku data
LEVELS = 176;
RES = 0.125;      % bin to bin resolution in km

%(Does not work for the whole globe! I.e., MinLongEast can never be equal to MaxLongEast.)
MinLongEast =  -85;  %  coordinates of whole area of interest: M2 (Cent Amer)
MaxLongEast =  -65;  %  
MINLATEAST  =   -5;  %
MAXLATEAST  =   20;  % 
%Do not set your southern border between -60 and -65. You can set the border at -65.

%With changing boxsize you can exclude information that produces boxes smaller than a
%specific area from being produced. The unit of boxsize is degree^2. Files
%whos surrounding lat long box is smaller or equal to boxsize are not being
%created. If you set boxsize to zero (default), all files are created whos horizontal
%projection has two or more contiguous pixels.
boxsize=0;

%%%outputdirectory
out_dir='/home/disk/bob/gpm/m2_cmed_ku/classify/ex_data/';

%%%inputdirectory for GPM 2Ku
main_dir='/home/disk/bob/gpm/';
in_dir=[main_dir,'extracted/2Ku/'];
coord_dir=[main_dir,'coord_Orbits_ku/'];

%%%temporal file folder for decompressing files
dirTEMP='/virtual/';

%%%flagEcho flag.  If set to 1, include flagEcho in processing.
useFlagEcho = 1;

%The 2Ku data should be the MODified version extracted and
%compressed (file.hdf.gz) from the original files
%and should be located in subfolders with years and months as substructure
%%The latent heating data should also be located in subfldersfiles, 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if MINLATEAST>-65 && MINLATEAST<-60
    error('Do not set the southern border between -60 and -65 but at -65.')
end

GRIDBOXSTEP = 0.05;
RAYS = 49;
if mod(RAYS,2) == 0
    error('RAYS must be an odd number in this code.')
end
SCAN_ANGLE = 17.0;   % +/-17 deg crosstrack
fileType = 'GPM-2Ku';
angleRANGE = (-SCAN_ANGLE:(SCAN_ANGLE/fix(RAYS/2)):SCAN_ANGLE); 
MinLongEastOrig=MinLongEast;
MaxLongEastOrig=MaxLongEast;

src_years=[2014];

[dummy, n_years]=size(src_years);
for yyy=1:n_years
    y_file=src_years(1,yyy);   

    %here I select the months
    if y_file==2014
        month_ini=05;
        month_fin=05;
    else
        month_ini=09;
        month_fin=09;
    end
    
    for mmm=month_ini:month_fin
        if mmm<10
            m_month=['0',num2str(mmm)];
        else
            m_month=[num2str(mmm)];
        end

        %%list of files in 2Ku and coordinates
        fileKu=dir([coord_dir,num2str(y_file),'/',m_month,'/*.hdf']);  
        files2Ku=char(fileKu(:).name);
        [n_files, lght]=size(files2Ku);
        %n_files=double(12); SRB
        
        %if y_file==2014 && mmm==03
        %    initial_file=1;
        %else
            initial_file=1;
        %end    

        for filenumber=initial_file:n_files 
            clear correctZ*
            clear geolocation*
            clear orbit*
            clear rainTypes*
            clear nearSurfRain*   
            %clear scanTime*
            clear dayFile*
            clear firstTime*
            clear lastTime*
            clear month*
            clear year*
  
            namein2Ku = char(files2Ku(filenumber,:));
            fprintf('\n')
            namein2Kunew=namein2Ku(1:length(namein2Ku)-10);
            disp(['n=',num2str(filenumber),'/',num2str(n_files),'  file analyze ',namein2Kunew]);
  
            orbit0 = str2double(namein2Ku(48:53));   %this is the orbit number     
            orbitS0 = namein2Ku(48:53);              %string of orbit number

            %asing variables with the string of dates
            dayFile0 = str2double(namein2Ku(29:30));
            month0 = str2double(namein2Ku(27:28));
            year0 = str2double(namein2Ku(23:26));
  
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Read the Coordinate Files to only extract data within range
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            filenameCoord=[coord_dir,num2str(y_file),'/',m_month,'/',namein2Kunew,'Coords.hdf'];
            sds_index = 1;   %(sdsID=1 latitude  in 2KuCoords)   sdsID starts in 0!
            latitude  = loadHDF_new(filenameCoord, sds_index);
            sds_index = 0;   %(sdsID=0 longitude in 2KuCoords)
            longitude = loadHDF_new(filenameCoord, sds_index);
            hdfml('closeall')
            
%figure(1); clf; 
%plot(longitude,latitude)

            losque=find(latitude>=MINLATEAST   & latitude<=MAXLATEAST & ...
                longitude>=MinLongEast & longitude<=MaxLongEast);
      
            [ll1, ll2]=size(losque);   %dimensions of array
            em=0;
            if ll1<2   % changed 12/01/2014 from ==0 to <2
                disp('Orbit is not in area.')
                em=1;  %this will tell us if there is no orbits within the area
            end
            clear losque latitude longitude sds_index ll1 ll2 filenameCoord
      
            if em==0   %em=1 will tell us if there is no orbits within the area
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % The following code loads the variables of interest.
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                gunzip([in_dir,num2str(y_file),'/',m_month,'/',namein2Kunew,'mod.hdf.gz'],dirTEMP);
                filename=[dirTEMP,namein2Kunew,'mod.hdf'];
 
                %% whos('geolocation0','rainTypes0','correctZ0','scanTimeAll0')
                %*****************************************************
                %here open the extracted files (mod files) by mzuluaga
                %*****************************************************
                sds_index = 0;   %(sdsID=0 longitude in 2Kumodified)  sdsID starts in 0!
                longitude = loadHDF_new(filename, sds_index);
                sds_index = 1;   %(sdsID=1 latitude in 2Kumodified)
                latitude  = loadHDF_new(filename, sds_index);
    
                [x1, x2]=size(latitude);
                geolocation0 = single(zeros(2,x1,x2));
                geolocation0(1,:,:) = latitude(:,:);
                geolocation0(2,:,:) = longitude(:,:);
                clear x1 latitude longitude

%figure(1); clf; 
%lons=reshape(geolocation0(2,:,:),1,[]);
%lats=reshape(geolocation0(1,:,:),1,[]);
%plot(lons,lats)
%hold on

                %% this is substituted instead of findcuts @@@@@@@@@@@@@@@@
                [x1, x2, x3]=size(geolocation0);   %dimensions of arrayexit
                wholeLongs0=zeros(x2,x3);
                wholeLongs0(:,:)=geolocation0(2,:,:);   
                wholeLats0=zeros(x2,x3);
                wholeLats0(:,:)=geolocation0(1,:,:);
                longvector=wholeLongs0;
                longvector(longvector<MinLongEast)=nan;
                longvector(longvector>MaxLongEast)=nan;
                longvector=max(longvector); %here a vector with longitudes of orbits in study area
                latvector=wholeLats0;
                latvector(latvector<MINLATEAST)=nan;
                latvector(latvector>MAXLATEAST)=nan;
                latvector=max(latvector);   %here a vector with latitudes of orbits in  study area
                longvector(isnan(latvector))=nan;

                sds_index = 2;   %(sdsID=2 correctZFactor in 2Kumodified)
                correctZ0 = loadHDF_new(filename,sds_index);
	
                sds_index = 3;   %(sdsID=3 nearSurfRain in 2Kumodified)
                nearSurfRain0 = loadHDF_new(filename,sds_index);
    
                sds_index = 4;   %(sdsID=4 nearSurfPhase in 2Kumodified)
                nearSurfPhase0 = loadHDF_new(filename,sds_index);
    
                sds_index = 8;   %(sdsID=8 hour in 2Kumodified)
                hour0 = loadHDF_new(filename,sds_index);
                hour0=single(hour0);
                sds_index = 9;   %(sdsID=9 minute in 2Kumodified)
                minute0 = loadHDF_new(filename,sds_index);
                minute0 = single(minute0);
                sds_index = 10;   %(sdsID=10 second in 2Kumodified)
                second0 = loadHDF_new(filename,sds_index);
                second0 = single(second0);
      
                [xx2, xx3]=size(hour0);
                firstTime0 = hour0(1)*3600.+minute0(1)*60+second0(1);
                lastTime0 = hour0(xx2)*3600+minute0(xx2)*60+second0(xx2);
	
                sds_index = 11;   %(sdsID=11 time in 2Kumodified)
                time0 = loadHDF_new(filename,sds_index);
                scanTimeAll0 = double(time0(4,:))*3600. + double(time0(5,:))*60. + double(time0(6,:));
                %scanTimeAll0 = scanTimeAll0';   %%%%'
    
                sds_index = 12;   %(sdsID=12 rainTypes in 2Kumodified)
                rainTypes0 = loadHDF_new(filename,sds_index);
	
                sds_index = 13;   %(sdsID=13 flagShallowRain in 2Kumodified)
                shallowRainTypes0 = loadHDF_new(filename,sds_index);
	
                sds_index = 14;   %(sdsID=14 widthBB in 2Kumodified)
                widthBB0 = loadHDF_new(filename,sds_index);
	
                sds_index = 15;   %(sdsID=15 heightBB in 2Kumodified)
                heightBB0 = loadHDF_new(filename,sds_index);
	
                if useFlagEcho
                    sds_index = 16;   %(sdsID=16 flagEcho in 2Kumodified)
                    flagEcho0 = loadHDF_new(filename,sds_index);
                end
                
                hdfml('closeall')
                delete(filename);    %%delete the decompressed files

                %%% *******************************************************

                %% whos('geolocation0','correctZ0','nearSurfRain0','nearSurfPhase0',)
                %% whos('scanTimeAll0','rainTypes0','shallowRainTypes0','widthBB0','heightBB0','flagEcho0')
                %%%here is selected only data on the study region
                geolocation0=geolocation0(:,:,isfinite(longvector));  %trim the geolocation
                correctZ0=correctZ0(:,:,isfinite(longvector));
                nearSurfRain0=nearSurfRain0(:,isfinite(longvector));  %%%%%%%%
                nearSurfPhase0=nearSurfPhase0(:,isfinite(longvector));  %%%%%%%%
                hour0=hour0(isfinite(longvector),1);
                minute0=minute0(isfinite(longvector),1);
                second0=second0(isfinite(longvector),1);
                scanTimeAll0=scanTimeAll0(:,isfinite(longvector));
                rainTypes0=rainTypes0(:,isfinite(longvector));
                shallowRainTypes0=shallowRainTypes0(:,isfinite(longvector));
                widthBB0=widthBB0(:,isfinite(longvector));
                heightBB0=heightBB0(:,isfinite(longvector));
                if useFlagEcho
                    flagEcho0=flagEcho0(:,:,isfinite(longvector));
                end
                
%lons=reshape(geolocation0(2,:,:),1,[]);
%lats=reshape(geolocation0(1,:,:),1,[]);
%plot(lons,lats,'r')

                clear filename namein2Ku sds_index
                %clear filename2H25 file2H25
                clear longvector latvector wholeLongs0 wholeLats0
                clear x1 x2 x3 xx2 xx3
	
                %Now chechk for holes in the orbit
                [x1, x2, x3]=size(geolocation0);
                wholeLongs0=zeros(x2,x3);
                wholeLongs0(:,:)=geolocation0(2,:,:);
                wholeLats0=zeros(x2,x3);
                wholeLats0(:,:)=geolocation0(1,:,:);
    
                long4=wholeLongs0(1,:);
                long5=zeros(1,length(long4)+1);
                long5(:,2:length(long5))=long4;
                long5(1,1)=long4(1,1);
                long4(1,length(long4)+1)=long4(1,length(long4));
	
                diffTest=abs(long4-long5);
                maximumTest=max(diffTest);
                
                % This is necessary if there are holes in orbit
                if maximumTest>0.15 
                    
                    %disp(['maximumTest = ',num2str(maximumTest),' (>0,15).  There are holes in the orbit.'])  % SRB
                    flag_hole=1;
                    donde = find(diffTest>=0.15);  % locates longitude range of the whole area of interest in the swath
                    orbit=num2str(orbit0);
                    orbitS=orbitS0;
                    dayFile=dayFile0;
                    month=month0;
                    year=year0;
                    geolocation=cat(3,geolocation0(:,:,1:donde-1));
                    correctZ=cat(3,correctZ0(:,:,donde:end),correctZ0(:,:,1:donde-1));
                    nearSurfRain=cat(2,nearSurfRain0(:,donde:end),nearSurfRain0(:,1:donde-1));
                    nearSurfPhase=cat(2,nearSurfPhase0(:,donde:end),nearSurfPhase0(:,1:donde-1));
                    scanTimeAll=cat(2,scanTimeAll0(1,donde:end),scanTimeAll0(1,1:donde-1));
                    firstTime=firstTime0;
                    lastTime=lastTime0;
                    rainTypes=cat(2,rainTypes0(:,donde:end),rainTypes0(:,1:donde-1));
                    shallowRainTypes=cat(2,shallowRainTypes0(:,donde:end),shallowRainTypes0(:,1:donde-1));
                    widthBB=cat(2,widthBB0(:,donde:end),widthBB0(:,1:donde-1));
                    heightBB=cat(2,heightBB0(:,donde:end),heightBB0(:,1:donde-1));
                    if useFlagEcho
                        flagEcho=cat(3,flagEcho0(:,:,donde:end),flagEcho0(:,:,1:donde-1));
                    end
                    
% $$$ %lons1=reshape(geolocation(2,:,:),1,[]);
% $$$ %lats1=reshape(geolocation(1,:,:),1,[]);
% $$$ %plot(lons1,lats1,'r')
% $$$ %fid55=fopen([out_dir,'test/orbits_with_holes.txt'],'a');
% $$$ %fprintf(fid55,'Hole in orbit %8.1f \n',orbit0);
% $$$ %fclose(fid55);

                else
                    
                    %disp(['maximumTest = ',num2str(maximumTest),' (<=0,15).  There are no holes in the orbit.'])  % SRB
                    flag_hole=0;
                    orbit=num2str(orbit0);
                    orbitS=orbitS0;
                    dayFile=dayFile0;
                    month=month0;
                    year=year0;
                    geolocation=geolocation0;
                    correctZ=correctZ0;
                    nearSurfRain=nearSurfRain0;
                    nearSurfPhase=nearSurfPhase0;
                    scanTimeAll=scanTimeAll0;
                    firstTime=firstTime0;
                    lastTime=lastTime0;
                    rainTypes=rainTypes0;
                    shallowRainTypes=shallowRainTypes0;
                    widthBB=widthBB0;
                    heightBB=heightBB0;
                    if useFlagEcho
                        flagEcho=flagEcho0;
                    end
                    
                end
                
                clear orbit0 orbitS0 dayFile0 month0 year0
                clear geolocation0 correctZ0 nearSurfRain0 nearSurfPhase0 scanTimeAll0
                clear firstTime0 lastTime0
                clear rainTypes0 shallowRainTypes0 widthBB0 heightBB0 
                if useFlagEcho
                    clear flagEcho0
                end
                clear directory
                clear namein2Ku
                clear sds_index
                %% whos('geolocation','correctZ','nearSurfRain','nearSurfPhase',)
                %% whos('scanTimeAll','rainTypes','shallowRainTypes','widthBB','heightBB','flagEcho')

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % The following code creates latitude & longitude vectors
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                [a1, a2, a3]=size(geolocation);
                SCANS = a3;
                clear a1 a2 a3
                clear latitudes longitudes
	
                latitudes(1:RAYS,1:SCANS)= geolocation(1,1:RAYS,1:SCANS);
                longitudes(1:RAYS,1:SCANS)= geolocation(2,1:RAYS,1:SCANS);
                
%%%%here it plot the orbit cutted 
%%plot(reshape(longitudes,1,[]),reshape(latitudes,1,[]),'m')
	
                clear geolocation SCANS
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                mm = longitudes >= MinLongEast & longitudes <= MaxLongEast;
                nn = find(sum(mm)>0);  % locates longitude range of the whole area of interest in the swath
                NN = length(nn);
                lat = double(latitudes(:,nn));
                long1 = double(longitudes(:,nn));
                
%plot(reshape(long1,1,[]),reshape(lat,1,[]),'g')
%hold on
                clear mm
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %The following code searches for data in the swath.
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                refl=correctZ(:,:,nn);
                reflpos=zeros(LEVELS,RAYS,NN);
                reflpos(refl>0)=1;    %assignate 1 where Zfact is > 0 avoid missing data
	
                clear nn NN refl
	
                reflpos(:,lat<MINLATEAST)=0;
                reflpos(:,lat>MAXLATEAST)=0;
	
                emptyscans=sum(sum(reflpos));  %plot(squeeze(emptyscans))
                outin=zeros(size(emptyscans)); %plot(squeeze(outin))
	
                if emptyscans==0   %here start to analyze scans if 0 jump to next orbit   
                    clear reflpos
                    clear emptyscans
                    clear long1 outin
                else
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %here assignates a variable called outin with type of data!!!
                    %search for data (beginning of contiguous scans with data get a 1 in
                    %outin, middle scans get a 2 and end scans a 3, non contiguous scans
                    %get a 4, empty scans are 0.
	  
                    if emptyscans(:,:,1)>0 && emptyscans(:,:,2)>0 %first scan
                        outin(1)=1;
                    elseif emptyscans(:,:,1)>0 && emptyscans(:,:,2)==0
                        outin(1)=4;
                    end
	  
                    if emptyscans(:,:,length(emptyscans)-1)>0 && emptyscans(:,:,length(emptyscans))>0 %last scan
                        outin(length(emptyscans))=3;
                    elseif emptyscans(:,:,length(emptyscans)-1)==0 && emptyscans(:,:,length(emptyscans))>0
                        outin(length(emptyscans))=4;
                    end
	  
                    for i=2:length(emptyscans)-1 %all scans in between first and last
                        if emptyscans(:,:,i)>0
                            if emptyscans(:,:,i-1)>0 && emptyscans(:,:,i+1)>0
                                outin(i)=2;
                            elseif emptyscans(:,:,i-1)>0 && emptyscans(:,:,i+1)==0
                                outin(i)=3;
                            elseif emptyscans(:,:,i-1)==0 && emptyscans(:,:,i+1)>0
                                outin(i)=1;
                            elseif emptyscans(:,:,i-1)==0 && emptyscans(:,:,i+1)==0
                                outin(i)=4;
                            end
                        end
                    end
                    
% $$$ figure(3); clf; 
% $$$ cos1=find(outin==1); [d1 d2]=size(cos1);
% $$$ plot(reshape(long1(:,cos1),1,RAYS*d1),reshape(lat(:,cos1),1,RAYS*d1),'g')
% $$$ hold on
% $$$ cos1=find(outin==2); [d1 d2]=size(cos1);
% $$$ plot(reshape(long1(:,cos1),1,RAYS*d1),reshape(lat(:,cos1),1,RAYS*d1),'r')
% $$$ cos1=find(outin==3); [d1 d2]=size(cos1);
% $$$ plot(reshape(long1(:,cos1),1,RAYS*d1),reshape(lat(:,cos1),1,RAYS*d1),'b')
% $$$ cos1=find(outin==4); [d1 d2]=size(cos1);
% $$$ plot(reshape(long1(:,cos1),1,RAYS*d1),reshape(lat(:,cos1),1,RAYS*d1),'black')
	
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %The following code calculates the lat and long boxes where there is
                    %data and cuts these boxes out of the swath.
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    reflpos2D=sum(reflpos);         %here accounts only for horizontal 2D       
                    contiguous=zeros(size(outin));  %var to account for single contigouos scans within swath    
	  
                    clear reflpos emptyscans
	  
                    indices=find(outin==4); %non contiguous scans with just a single full pixel are ignored
                    [e1, e2]=size(indices);
                    for i=1:e1
                        a=indices(i);
                        for j=1:RAYS-1
                            if reflpos2D(:,j,a)>0 && reflpos2D(:,j+1,a)>0
                                contiguous(:,:,a)=1;
                            end
                        end
                    end
	  
                    clear reflpos2D indices e1 e2 a
	  
                    %Localize the outerbounds of the regions where data is located
                    %search for the longitude limits
                    longminimum=zeros(size(outin));
                    longminimum(:,:,:)=nan;
                    longminimum(1,1,outin==1)=min(long1(:,outin==1));
                    longminimum(1,1,contiguous==1)=min(long1(:,contiguous==1));
                    longminimum=longminimum(isfinite(longminimum));
	
                    longmaximum=zeros(size(outin));
                    longmaximum(:,:,:)=nan;
                    longmaximum(1,1,outin==3)=max(long1(:,outin==3));
                    longmaximum(1,1,contiguous==1)=max(long1(:,contiguous==1));
                    longmaximum=longmaximum(isfinite(longmaximum));
	  
                    longlimits2=cat(1,longminimum,longmaximum);
                    [q1, q2, q3]=size(longlimits2);
                    clear longlimits
                    longlimits=zeros(2,q3);
                    longlimits(:,:)=longlimits2;
	  
                    clear long1 outin contiguous
                    clear longminimum longmaximum longlimits2
	
                    %put small boxes together
                    if q3>1   %so What if q3 is less than 1????
                        for i=1:q3
                        %for i=7:7   SRB
                            if longlimits(1,i)>=longlimits(2,i)
                                [d1, d2]=size(longlimits);
                                fprintf('\n')
                                disp('There is a reverse maximum in swath')
                                while longlimits(2,i+1)<longlimits(1,i)
                                    longlimits(1,i)=longlimits(1,i+1);
                                    longlimits(:,i+1)=[];
                                    [d1, d2]=size(longlimits);
                                    if i+1>d2
                                        break
                                    end
                                end
                            else
                                [d1, d2]=size(longlimits);
                                while longlimits(1,i+1)<longlimits(2,i)
                                    longlimits(2,i)=longlimits(2,i+1);
                                    longlimits(:,i+1)=[];
                                    [d1, d2]=size(longlimits);
                                    if i+1>d2
                                        break
                                    end
                                end
                            end
                            if i+2>d2
                                break
                            end
                        end
                    end
                    clear q1 q2 q3 d1 d2
	  
                    %longlimits has now the longitude limites of the boxes.
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %The following code goes now through all boxes listed in longlimits and
                    %calculates all needed output.
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    clear f1 f2
                    [f1, f2]=size(longlimits);
	  
                    fprintf('\n')
                    disp(['Number of boxes in swath to be analyzed: ',num2str(f2)])
	  
                    if f2>0        %for more than 1 region!!!
                        for inoneorbit=1:f2
                            MINLONGEAST=longlimits(1,inoneorbit);
                            MAXLONGEAST=longlimits(2,inoneorbit);
	      
                            if MINLONGEAST<MinLongEast
                                MINLONGEAST=MinLongEast;
                            end
	      
                            if MAXLONGEAST>MaxLongEast
                                MAXLONGEAST=MaxLongEast;
                            end
	      
                            if MAXLONGEAST<(MinLongEast+0.05) || MINLONGEAST>(MaxLongEast-0.05)
                                fprintf('\n')
                                disp('box found in a border of the region...skiping.this.box...')
                                %fprintf('\n')
                                clear MINLONGEAST
                                clear MAXLONGEAST
                            else
		
                                m = longitudes >= MINLONGEAST & longitudes <= MAXLONGEAST;
                                n = find(sum(m)>0);  % locates appropriate longitude range
                                N = length(n);
                                N1 = n(1,1);  % index of the first good value in scans
                                Nn = n(1,N);  % index of the last good value in scans
% $$$ figure(4); clf;
% $$$ cosM=find(m==1);
% $$$ plot(longitudes(cosM),latitudes(cosM))
% $$$ hold on
% $$$ lons2=longitudes(:,n); lats2=latitudes(:,n);
% $$$ plot(reshape(lons2,1,RAYS*N),reshape(lats2,1,RAYS*N),'r')
		
                                clear m
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                centreScan = ceil((Nn-N1)/2);
                                if Nn==N1
                                    centreScan=1;
                                end
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                %here select lats,lons and rain type within the analized box
                                lat = double(latitudes(:,n));   
                                long = double(longitudes(:,n));
                                sR1 = double(nearSurfRain(:,n));     %nearSurfRain
                                sP1 = double(nearSurfPhase(:,n));    %nearSurfPhase
                                rt1 = double(rainTypes(:,n));        %rainTypes
                                srt1 = double(shallowRainTypes(:,n));%shallowRainTypes
                                wBB1 = double(widthBB(:,n));         %widthBB
                                hBB1 = double(heightBB(:,n));        %heightBB
                                if useFlagEcho
                                    fEcho1 = double(flagEcho(:,:,n));    %flagEcho
                                end
                                
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                %here is selected the scantime of the box 
                                scanTime = scanTimeAll(N1:Nn); 
                                %Calculates the correct time and date (UTC)
                                groundTrackDate=dateSorterGPM(centreScan,dayFile,month,year,firstTime,lastTime,scanTime);
		
                                %function to calculate the date and time in particular  format
                                date_box=datestring(groundTrackDate,1);
		
                                clear centreScan scanTime
		
                                %here it assignates the file name of the analyzed box
                                fileInfo = [fileType,'.',orbitS,'.',date_box,'.'];
                                clear groundTrackDate
		
                                fprintf('\n')
                                disp(['the analyzed box will be named ',fileInfo])
		
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                %The following code creates an array of the radar reflectivity at all height bins (Z##)
                                %It does no create array, it creates LEVELS variables with Z at each level
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		
                                for k=1:LEVELS;  % LEVELS = 176 for REAL files
                                    eval(['Z',num2str(k),'(1:RAYS,1:N)=correctZ(' num2str(k),',1:RAYS,N1:Nn);']);
                                    eval(['Z',num2str(k),' = double(Z' num2str(k),');']);
                                end

                                clear N1 Nn
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                % The following code corrects the geolocation and gives the latitude & longitude
                                % coordinates based on scan & height.  NOTE: Level 176 is the earth ellipsoid
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                heightmatrix=zeros(RAYS,N,LEVELS);    %This is for Reflectivity
                                for i=1:LEVELS
                                    heightmatrix(:,:,i)=(LEVELS-i)*RES; %height in km   from 22km to 0
                                end

                                angle2D=zeros(RAYS,N);
                                for i=1:N
                                    angle2D(:,i)=angleRANGE;   %here assignate the angle of each scan in 2D    
                                end
		
                                anglerangematrix=zeros(size(heightmatrix));
                                for i=1:LEVELS
                                    anglerangematrix(:,:,i)=angle2D;
                                end

                                %matrix with shift distance for each  level       
                                kmshift=heightmatrix.*tan(deg2rad(anglerangematrix));  %%in km 
                                range=abs(km2deg(kmshift));                            %% in deg  

                                %calculate the angle in which the geolocation needs to be shifted
                                az = azimuth(lat(1,1:N),long(1,1:N),lat(RAYS,1:N),long(RAYS,1:N));
                                directionshiftBOTTOM=az;
                                if az < 180
                                    directionshiftTOP = az + 180;
                                else
                                    directionshiftTOP = az - 180;
                                end
		
                                lat3D=zeros(size(heightmatrix));   %this is latitudes in 3D matrix   
                                for i=1:LEVELS
                                    lat3D(:,:,i)=lat;
                                end

                                long3D=zeros(size(heightmatrix));  %this is longitudes in 3D matrix 
                                for i=1:LEVELS
                                    long3D(:,:,i)=long;
                                end

                                clear heightmatrix angle2D anglerangematrix kmshift az

                                %this is a matrix with 24 values of N-scans azimuths      
                                directionshiftBOTTOM2D=zeros(fix(RAYS/2),N);
                                for i=1:fix(RAYS/2)
                                    directionshiftBOTTOM2D(i,:)=directionshiftBOTTOM;
                                end
		
                                directionshiftTOP2D=zeros(fix(RAYS/2+1),N);
                                for i=1:fix(RAYS/2+1)
                                    directionshiftTOP2D(i,:)=directionshiftTOP;
                                end
		
                                directionshiftBOTTOM3D=zeros(fix(RAYS/2),N,LEVELS);
                                for i=1:LEVELS
                                    directionshiftBOTTOM3D(:,:,i)=directionshiftBOTTOM2D;
                                end

                                directionshiftTOP3D=zeros(fix(RAYS/2+1),N,LEVELS);
                                for i=1:LEVELS
                                    directionshiftTOP3D(:,:,i)=directionshiftTOP2D;
                                end

                                %%%% lats and lons corrected for Reflectivity
                                [latsB,longsB] = reckon(lat3D(1:fix(RAYS/2),1:N,1:LEVELS),long3D(1:fix(RAYS/2),1:N,1:LEVELS),...
                                    range(1:fix(RAYS/2),1:N,1:LEVELS),directionshiftBOTTOM3D(:,:,1:LEVELS));
		
                                [latsT,longsT] = reckon(lat3D(fix(RAYS/2+1):RAYS,1:N,1:LEVELS),long3D(fix(RAYS/2+1):RAYS,1:N,1:LEVELS),...
                                    range(fix(RAYS/2+1):RAYS,1:N,1:LEVELS),directionshiftTOP3D(:,:,1:LEVELS));

                                clear range directionshift* lat3D long3D
		
                                lats=cat(1,latsB,latsT);     %corrected latitude
                                longs=cat(1,longsB,longsT);  %corrected longitude

%% plot reflectivity at every level (earth ellipsoid at Z176)
%  name = 'gpm.20140705.205454.';
%  name = 'gpm.20140706.200809.';
%  for z=169:176
%      zstr=num2str(z);
%      km = (176 - z) * 0.125;
%      eval(['figure(',zstr,')']); clf;
%      eval(['R',zstr,' = Z',zstr,';']);
%      eval(['ind = find(R',zstr,'==-9999);']);
%      eval(['R',zstr,'(ind)=NaN;']);
%      eval(['R',zstr,' = 10*log10(R',zstr,');']);
%      eval(['contourf(longs(:,:,',zstr,'),lats(:,:,',zstr,'),R',zstr,')']);
%      colorbar
%      %caxis([30,38])
%      xlabel('longitude')
%      ylabel('latitude')
%      title(['Reflectivity at ',num2str(km),' km'])
%      %axis([-81,-77,29.5,34.5])
%      %axis([-56,-50,-42,-35])   % 20140705
%      axis([-43,-35,-48,-42])   % 20140706
%      colormap jet
%      %eval(['minR = nanmin(nanmin(R',zstr,'));']);
%      eval(['maxR = nanmax(nanmax(R',zstr,'));']);
%      %str={['min = ',num2str(minR],['max = ',num2str(maxR)]};
%      str={['max = ',num2str(maxR)]};
%      %annotation('textbox',[0.15,0.75,0.1,0.1],'String',str);
%      annotation('textbox',[0.5,0.75,0.1,0.1],'String',str);
%      if length(zstr)==1
%          zstr=['00',zstr];
%      elseif length(zstr)==2
%          zstr=['0',zstr];
%      end
%      eval(['print -dpng ',name,zstr,'.png']);
%  end

                                clear latsB latsT longsB longsT
                                
%% whos('lats','longs','lats_LH','longs_LH')	      
% $$$ figure(5); clf;   %plot the corrected coordinates
% $$$ plot(long(:,1),lat(:,1),'--rs')                   %this are the original coordinates
% $$$ hold on
% $$$ plot([long(1,1),long(RAYS,1)],[lat(1,1),lat(RAYS,1)],'-ks') %straigt line between edges   
% $$$ plot(longs(:,1,1),lats(:,1,1),'--gs')           %this are the corrected 21.875km coord
% $$$ plot(longs(:,1,LEVELS),lats(:,1,LEVELS),'--s')         %this are the corrected 0.00km coord

% $$$ plot(long(:,1),lat(:,1),'--rs')                   %this are the original coordinates
% $$$ plot([long(1,1),long(RAYS,1)],[lat(1,1),lat(RAYS,1)],'-ks') %straigt line between edges   
% $$$ plot(longs_LH(:,1,1),lats_LH(:,1,1),'--gs')           %this are the corrected 0.00km coord
% $$$ plot(longs_LH(:,1,19),lats_LH(:,1,19),'--s')         %this are the corrected 17km coord


                                % The product of this section of code
                                % is the corrected latitude and longitude called lats## and longs##
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                % The following code calculates the dimensions of the box in
                                % the new cartesian grid on which the reflectivity will be
                                % interpolated. Gridspacing is 0.05 degrees.
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		
                                reflsmall=correctZ(:,:,n);    %select Zval within the analyzed box   
                                reflpossmall=zeros(LEVELS,RAYS,N);  %create a matrix to find where Zval data exist          
                                reflpossmall(reflsmall>0)=1;
                                reflpos2Dsmall=sum(reflpossmall); %where in 2Dvertical dimension data are     
		
                                clear reflsmall reflpossmall
		
                                minlatND = min(min(lat(reflpos2Dsmall>0)));  %to find the latitude limits    
                                maxlatND = max(max(lat(reflpos2Dsmall>0)));         
		
                                modMinLat = mod(minlatND, 0.05);         
                                modMaxLat = mod(maxlatND, 0.05);         
            
                                minlat = minlatND - modMinLat;         
                                maxlat = maxlatND + (0.05 - modMaxLat);         
		
                                clear minlatND maxlatND modMinLat modMaxLat
		
                                if minlat < MINLATEAST
                                    minlat = MINLATEAST;
                                end
		
                                if maxlat > MAXLATEAST
                                    maxlat = MAXLATEAST;
                                end
		
                                minlongND = min(min(long(reflpos2Dsmall>0)));  %to find the longitude limits
                                maxlongND = max(max(long(reflpos2Dsmall>0)));
		
                                clear reflpos2Dsmall
		
                                modMinLong = mod(minlongND, 0.05);
                                modMaxLong = mod(maxlongND, 0.05);
		
                                modMINLONGEAST = mod(MINLONGEAST, 0.05);
                                modMAXLONGEAST = mod(MAXLONGEAST, 0.05);
		
                                minlong = minlongND - modMinLong;
                                maxlong = maxlongND + (0.05 - modMaxLong);
		
                                clear minlongND maxlongND modMinLong modMaxLong
		
                                if minlong<MINLONGEAST-modMINLONGEAST
                                    minlong=MINLONGEAST-modMINLONGEAST;
                                end
		
                                if maxlong>MAXLONGEAST+(0.05 - modMAXLONGEAST)
                                    maxlong=MAXLONGEAST+(0.05 - modMAXLONGEAST);
                                end
		
                                if minlong < MinLongEast
                                    minlong = MinLongEast;
                                end
		
                                if maxlong > MaxLongEast
                                    maxlong = MaxLongEast;
                                end
		
                                clear modMINLONGEAST modMAXLONGEAST MINLONGEAST MAXLONGEAST
		
                                boundaries=[minlong maxlong minlat maxlat];
                                boxsize1=(maxlong-minlong)*(maxlat-minlat);
		
                                %fprintf('\n')
                                disp(['Box: ',num2str(boundaries)])
		
                                clear boundaries
		
                                [longitudeGridEast,latitudeGridEast] = ...
                                    meshgrid(minlong:GRIDBOXSTEP:maxlong,minlat:GRIDBOXSTEP:maxlat);

                                %%% whos('longitudeGridEast','latitudeGridEast')
                                %This is the grid in which the data will be interpolated (new grid).
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                %here it will create the files with the information of pixels with
                                %Z>0 systems found in the boxes previously analyzed      
                                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                if minlong<maxlong && minlat<maxlat && boxsize1>boxsize %check for boxsize at beggining
                                    clear boxsize1
                                    sizeGrid = size(longitudeGridEast);
                                    gridRows = sizeGrid(1,1);
                                    gridColumns = sizeGrid(1,2);
                                    numSwaths = N;
		  
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %This code creates a grid with the size of the new grid with 0 inside
                                    %the swath and 1 outside the swath. This is usefull because we     
                                    %dont have to interpolate outside the swath.
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    empty=zeros(size(longitudeGridEast));
		  
                                    for i=1:sizeGrid(1,1)
                                        for j=1:sizeGrid(1,2)
                                            temp=find(lat>latitudeGridEast(i,j)-0.1 & lat<latitudeGridEast(i,j)+0.1...
                                                & long>longitudeGridEast(i,j)-0.1 & long<longitudeGridEast(i,j)+0.1);
                                            if isempty(temp)
                                                empty(i,j)=1;
                                            end
                                            temp=[];
                                        end
                                    end
		  
                                    clear lat long sizeGrid temp  %%delete the non corrected lat and long
		  
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    % The following code creates column versions of the regularly 
                                    % gridded latitude & logitude coordinates.
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    gridLongitude = reshape(longitudeGridEast, numel(longitudeGridEast), 1);
                                    gridLatitude = reshape(latitudeGridEast, numel(latitudeGridEast), 1);
		  
                                    %%whos('gridLongitude','longitudeGridEast','gridLatitude','latitudeGridEast')
                                    clear longitudeGridEast latitudeGridEast
                                    
                                    noSwath = reshape(empty, numel(empty), 1);
                                    clear empty
		  
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    % The following code creates column versions of the reflectivity data at all levels.
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    sRLEVELSc = reshape(sR1, numel(sR1), 1);   %near surface rain for all scans%%%%%%%%
                                    sPLEVELSc = reshape(sP1, numel(sP1), 1);   %near surface phase for all scans%%%%%%%%
                                    rtLEVELSc = reshape(rt1, numel(rt1), 1);   %rain type for all scans         
                                    srtLEVELSc = reshape(srt1, numel(srt1), 1);%shallow rain type for all scans         
                                    wBBLEVELSc = reshape(wBB1, numel(wBB1), 1);%width bright band for all scans         
                                    hBBLEVELSc = reshape(hBB1, numel(hBB1), 1);%height bright band for all scans         
                                    if useFlagEcho
                                        fEchoLEVELSc = reshape(fEcho1, numel(fEcho1), 1);%flagEcho for all scans         
                                        clear fEcho1
                                    end
                                    clear sR1 sP1 rt1 srt1 wBB1 hBB1
                                        
                                    for level = 1:LEVELS
                                        eval(['Z',num2str(level),'c=reshape(Z',num2str(level),...
                                            ',prod(size(Z',num2str(level),')), 1);']);
                                    end
                                        
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    % Gets the surface rain & phase, the
                                    % rain type & shallow rain type, and
                                    % the width & height of the bright band
                                    % for the new grid by searching for the 
                                    % nearest neighbour in the old grid.
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    [d1, d2]=size(gridLongitude);
                                    gridSurfRain=zeros(d1,1);
                                    gridSurfPhase=zeros(d1,1);
                                    gridRainType=zeros(d1,1);
                                    gridShallowRainType=zeros(d1,1);
                                    gridWidthBB=zeros(d1,1);
                                    gridHeightBB=zeros(d1,1);
                                    for k = 1:d1;
                                        gridSurfRain(k,1) = getNearSurfRain(gridLongitude(k), gridLatitude(k),...
                                            reshape(longs(:,:,LEVELS),[],1), reshape(lats(:,:,LEVELS),[],1), sRLEVELSc); %%%%%%%%
                                        gridSurfPhase(k,1) = get2Dinterp(gridLongitude(k), gridLatitude(k),...
                                            reshape(longs(:,:,LEVELS),[],1), reshape(lats(:,:,LEVELS),[],1), sPLEVELSc);
                                        gridRainType(k,1) = getRaintype(gridLongitude(k), gridLatitude(k),...
                                            reshape(longs(:,:,LEVELS),[],1), reshape(lats(:,:,LEVELS),[],1), rtLEVELSc);
                                        gridShallowRainType(k,1) = get2Dinterp(gridLongitude(k), gridLatitude(k),...
                                            reshape(longs(:,:,LEVELS),[],1), reshape(lats(:,:,LEVELS),[],1), srtLEVELSc);
                                        gridWidthBB(k,1) = get2Dinterp(gridLongitude(k), gridLatitude(k),...
                                            reshape(longs(:,:,LEVELS),[],1), reshape(lats(:,:,LEVELS),[],1), wBBLEVELSc);
                                        gridHeightBB(k,1) = get2Dinterp(gridLongitude(k), gridLatitude(k),...
                                            reshape(longs(:,:,LEVELS),[],1), reshape(lats(:,:,LEVELS),[],1), hBBLEVELSc);
                                    end  %%%whos('rtLEVELSc','gridRainType','gridLongitude','gridLatitude')

                                    if useFlagEcho
                                        gridFlagEcho=zeros(d1*LEVELS,1);
                                        for kk = 1:LEVELS
                                            if rem(kk,10) == 0
                                                disp(['Interp flagEcho level: ',num2str(kk)])
                                            end
                                            for k = 1:d1
                                                gridFlagEcho(((kk-1)*LEVELS)+k,1) = get2Dinterp(gridLongitude(k), gridLatitude(k),...
                                                reshape(longs(:,:,kk),[],1), reshape(lats(:,:,kk),[],1), fEchoLEVELSc);
                                            end
                                        end
                                        clear fEchoLEVELSc
                                    end

                                    clear sRLEVELSc sPLEVELSc rtLEVELSc srtLEVELSc wBBLEVELSc hBBLEVELSc d1 d2

                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    % 1. Assign gridRainType values of 1 (stra),2 (conv), 3 (other)
                                    % 2. Add a new variable called gridRainTypeRaw which retain the orig gridRainType values
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    gridRainTypeRaw=gridRainType;
		  
                                    gridRainType(gridRainType>0 & fix(gridRainType/10000000.)==1) = 1;
                                    gridRainType(gridRainType>0 & fix(gridRainType/10000000.)==2) = 2;
                                    gridRainType(gridRainType>0 & fix(gridRainType/10000000.)==3) = 3;
                                        
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    % 1. Assign gridSurfPhase values of 1(solid),2 (mixed phase), 3 (liquid) to all values
                                    %    that are not negative or missingvals (255)
                                    % 2. Add a new variable called gridSurfPhaseRaw which retain the orig gridSurfPhase vals
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    gridSurfPhaseRaw=gridSurfPhase;
		  
                                    gridSurfPhase(gridSurfPhase>0 & gridSurfPhase~=255 & fix(gridSurfPhase/100)==0) = 0;
                                    gridSurfPhase(gridSurfPhase>0 & gridSurfPhase~=255 & fix(gridSurfPhase/100)==1) = 1;
                                    gridSurfPhase(gridSurfPhase>0 & gridSurfPhase~=255 & fix(gridSurfPhase/100)==2) = 2;
                                        
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %make reflectivity data ready for output file
                                    values = cat(1, Z176c, Z175c, Z174c, Z173c, Z172c, Z171c,...
                                        Z170c, Z169c, Z168c, Z167c, Z166c, Z165c, Z164c, Z163c, Z162c, Z161c,...
                                        Z160c, Z159c, Z158c, Z157c, Z156c, Z155c, Z154c, Z153c, Z152c, Z151c,...
                                        Z150c, Z149c, Z148c, Z147c, Z146c, Z145c, Z144c, Z143c, Z142c, Z141c,...
                                        Z140c, Z139c, Z138c, Z137c, Z136c, Z135c, Z134c, Z133c, Z132c, Z131c,...
                                        Z130c, Z129c, Z128c, Z127c, Z126c, Z125c, Z124c, Z123c, Z122c, Z121c,...
                                        Z120c, Z119c, Z118c, Z117c, Z116c, Z115c, Z114c, Z113c, Z112c, Z111c,...
                                        Z110c, Z109c, Z108c, Z107c, Z106c, Z105c, Z104c, Z103c, Z102c, Z101c,...
                                        Z100c, Z99c, Z98c, Z97c, Z96c, Z95c, Z94c, Z93c, Z92c, Z91c,...
                                        Z90c, Z89c, Z88c, Z87c, Z86c, Z85c, Z84c, Z83c, Z82c, Z81c,...
                                        Z80c, Z79c, Z78c, Z77c, Z76c, Z75c, Z74c, Z73c, Z72c, Z71c,...
                                        Z70c, Z69c, Z68c, Z67c, Z66c, Z65c, Z64c, Z63c, Z62c, Z61c,...
                                        Z60c, Z59c, Z58c, Z57c, Z56c, Z55c, Z54c, Z53c, Z52c, Z51c,...
                                        Z50c, Z49c, Z48c, Z47c, Z46c, Z45c, Z44c, Z43c, Z42c, Z41c,...
                                        Z40c, Z39c, Z38c, Z37c, Z36c, Z35c, Z34c, Z33c, Z32c, Z31c,...
                                        Z30c, Z29c, Z28c, Z27c, Z26c, Z25c, Z24c, Z23c, Z22c, Z21c,...
                                        Z20c, Z19c, Z18c, Z17c, Z16c, Z15c, Z14c, Z13c, Z12c, Z11c,...
                                        Z10c, Z9c,  Z8c,  Z7c,  Z6c,  Z5c,  Z4c,  Z3c,  Z2c,  Z1c);
		  
                                    % CHECK TO SEE IF THIS IS STILL VALID %
                                    badvalues = find(values < 0);
		  
                                    reflectivity = 10.^(values./1000); %% Z factor to reflectivity [dBZ -> mm6/m3]
                                    % this is 1000 because Z-fact is scaled  by 100.
                                    %and do the same for lats and longs
                                    fliplat=flipdim(lats,3);
                                    latitude=reshape(fliplat(:,:,1:LEVELS),[],1);
		  
                                    fliplong=flipdim(longs,3);
                                    longitude=reshape(fliplong(:,:,1:LEVELS),[],1);
		  
                                    reflectivity(badvalues) = -999;

                                    clear values badvalues fliplat fliplong
                                        
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    % saves the 19 output files
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    dir_out=[out_dir,date_box(1:4),'/',date_box(5:6),'/'];
                                    %dir_out=[out_dir,date_box(1:4),'/',date_box(5:6),'/TEST/'];
                                    if ~exist(dir_out,'dir')
                                        mkdir(dir_out)
                                    end
                                    
                                    fid = fopen([dir_out, fileInfo, 'hdfFileName.txt'],'w');
                                    fprintf(fid,'%s', namein2Kunew);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'numSwaths.txt'],'w');
                                    fprintf(fid,'%d', numSwaths);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'noSwath.txt'],'w');
                                    fprintf(fid,'%3.1f\t', noSwath);
                                    fclose(fid);

                                    fid = fopen([dir_out, fileInfo, 'longitude.txt'],'w');
                                    fprintf(fid,'%10.5f\t', longitude);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'latitude.txt'],'w');
                                    fprintf(fid,'%9.5f\t', latitude);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'reflectivity.txt'],'w');
                                    fprintf(fid,'%8.3f\t', reflectivity);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'gridRows.txt'],'w');
                                    fprintf(fid,'%d', gridRows);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'gridColumns.txt'],'w');
                                    fprintf(fid,'%d', gridColumns);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'gridLongitude.txt'],'w');
                                    fprintf(fid,'%10.5f\t', gridLongitude);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'gridLatitude.txt'],'w');
                                    fprintf(fid,'%9.5f\t', gridLatitude);
                                    fclose(fid);
		  
                                    fid = fopen([dir_out, fileInfo, 'gridNearSurfRain.txt'],'w');
                                    fprintf(fid,'%9.5f\t', gridSurfRain);
                                    fclose(fid);

                                    fid = fopen([dir_out, fileInfo, 'gridNearSurfPhase.txt'],'w');
                                    fprintf(fid,'%2.0f\t', gridSurfPhase);
                                    fclose(fid);

                                    fid = fopen([dir_out, fileInfo, 'gridNearSurfPhaseRaw.txt'],'w');
                                    fprintf(fid,'%4.0f\t', gridSurfPhaseRaw);
                                    fclose(fid);

                                    fid = fopen([dir_out, fileInfo, 'gridRainType.txt'],'w');
                                    fprintf(fid,'%6.0f\t', gridRainType);
                                    fclose(fid);

                                    fid = fopen([dir_out, fileInfo, 'gridRainTypeRaw.txt'],'w');
                                    fprintf(fid,'%9.0f\t', gridRainTypeRaw);
                                    fclose(fid);

                                    fid = fopen([dir_out, fileInfo, 'gridShallowRainType.txt'],'w');
                                    fprintf(fid,'%6.0f\t', gridShallowRainType);
                                    fclose(fid);

                                    fid = fopen([dir_out, fileInfo, 'gridWidthBB.txt'],'w');
                                    fprintf(fid,'%8.3f\t', gridWidthBB);
                                    fclose(fid);

                                    fid = fopen([dir_out, fileInfo, 'gridHeightBB.txt'],'w');
                                    fprintf(fid,'%8.3f\t', gridHeightBB);
                                    fclose(fid);
                                    
                                    if useFlagEcho
                                        fid = fopen([dir_out, fileInfo, 'gridFlagEcho.txt'],'w');
                                        fprintf(fid,'%8.3f\t', gridFlagEcho);
                                        fclose(fid);
                                    end

                                    clear date_box dir_out
                                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                else
                                    fprintf('\n')
                                    disp('Empty box or box too small.')
                                    clear boxsize1
                                end % end of boxsize restriction
                                
                                clear minlat
                                clear minlong
                                clear maxlat
                                clear maxlong
                                clear fid*
                                clear fileInfo
                                clear gridRows
                                clear gridColumns
                                clear gridLongitude
                                clear gridLatitude
                                clear gridRainType
                                clear gridShallowRainType
                                clear gridRainTypeRaw
                                clear gridSurfRain 
                                clear gridSurfPhase 
                                clear gridWidthBB
                                clear gridHeightBB
                                if useFlagEcho
                                    clear gridFlagEcho
                                end
                                clear i j k
                                clear lats
                                clear latitude
                                clear longs
                                clear longitude
                                clear n
                                clear noSwath
                                clear N
                                clear numSwaths
                                clear reflectivity
                                clear Z*
                                clear fx1 fx2
                            end % endif of restriction of limits being to close to the boundaries    
                        end  %endfor of analyzed regions
                    end   %endif of regions found having contigous data restriction for more than 1 region!!!
                end  %endif for orbits without correctZfact data in area of interest
                clear correctZ rainTypes shallowRainTypes
                clear nearSurfRain nearSurfPhase
                clear widthBB heightBB 
                if useFlagEcho
                    clear flagEcho
                end
                clear hour second minute scanTimeAll 
                clear latitudes longitudes
                clear namein2Kunew
            else
                fprintf('\n')
                disp(['The file ',namein2Kunew,' is not included in the limits']);
	
                clear filename namein2Ku* sds_index
                clear geolocation0 longvector latvector wholeLongs0 wholeLats0
                clear x1 x2 x3 em
            end    %endif for em=1 so no orbits within the area
        end  %end of for for different orbit files
    end  %loop through the different months
end  %loop through the different years

t=clock;
datestamper(t);
