function [] = netcdfwrite(file)
% 
% This function writes concatenated file to netcdf
% May, 2022
%
% Modified to include translation velocity, August, 2020
% Modified to include polyfilter output and scalar control paramters,
% November, 2020
% Revised May, 2022 to use with Octave and FAST
%-------------------------------------------------------------------------
%
load(file)
[n,m]=size(vstore);
%
[~,jmax]=min(vstore,[],2);
jmax=jmax-1;
dum=zeros(n,400);
yeartemp=zeros(1,jmax(1));
%
[~,jmax]=min(vstore,[],2);
jmax=jmax-1;
%
%  Calculate translation velocity and wind speed increments
%
if exist('u850store','var')
    [ut,vt,uinc,vinc ] = utransfull(latstore,longstore,vstore,u850store,v850store);
    if exist('vsestore','var')
        [~,~,uincse,vincse ] = utransfull(latstore,longstore,vsestore,u850store,v850store);
    end    
else
    [ut,vt,uinc,vinc ] = utransfull(latstore,longstore);
    if exist('vsestore','var')
        uincse=uinc;
        vincse=vinc;
    end
end
% 
vnet=vstore+sqrt(uinc.^2+vinc.^2);
[vnetmax,ivmax]=max(vnet,[],2);    % Storm lifetime maximum wind
vnetmax=transpose(vnetmax);
%
% Create mask for secondary eyewall events and find maximum ground relative
% wind for them
%
if exist('vsestore','var')
    semask=min(vsestore,1);
    vnetse=semask.*(vsestore+sqrt(uincse.^2+vincse.^2));
end 
if exist('T600store','var')
    [router] = roquick(T600store,rmstore,vstore);
end  
%
passi=zeros(n,m)+1;
if strcmp(shape,'poly')
    [polytype, passi, passe, xint, yint, jint, kint, kfrac] = polyfilter(latstore,longstore);
    if strcmp(polytype,'closed')
        poly=load('poly.in'); 
        x=poly(:,1); 
        y=poly(:,2);
        A=polyarea(x',y');
        if A < poly_area
            polytype='open';
        end
    end        
    %
    % Find maximum ground-relative wind at time track crosses line segements
    %
    [~,xsize]=size(xint);
    vseg=zeros(n,xsize); % Added 11/28/17
    for nx=1:xsize % Find intersection with maximum wind
        for i=1:n
            vseg(i,nx)=passe(i).*vnet(i,max(kint(i,nx),1)); % Changed 11/28/17
        end
    end   
    [vmax,iimax]=max(vseg,[],2);
    kmax2=zeros(1,n);
    for i=1:n
        kmax2(i)=kint(i,iimax(i));
    end  
    kmax2=max(kmax2,1);
    vmax=vmax';
    if strcmp(polytype,'closed')
        vmaxcl=max(vnet.*passi,[],2);  % Maximum ground-relative wind speed within closed polygon
        vmaxcl=transpose(vmaxcl);
        vmax=max(vmaxcl,vmax);         % For narrow closed polygons, max wind may be on boundary
        clear vmaxcl
    end      
else  % Case of circular filter, find maximum wind at point of interest
    if city_radius < 2000 % Exclude cases where city radius encompasses whole basin
        magfac=1.0;
        [~, vmax, kmax2] = best_filter(vnet,latstore,longstore,city_radius,clat,clong);
        [vs, dir1, dayk] = pointseriesn(latstore,longstore,vstore,rmstore.*magfac, ...
                 vsestore,rmsestore.*magfac,yearstore,monthstore,daystore,hourstore,uinc,vinc, ...
                 clat,clong,timeres);     
        vsmax=max(vs,[],2);
        vsmax=transpose(vsmax);
    else
        vmax=vnetmax;
        kmax2=ivmax;
    end    
end
[mf,nf]=size(vmax);
if mf > nf
    vmax=transpose(vmax);
end 
eset=file;
%
        model=char(model);
        gmeth=char(gmeth);
        vcrit=str2double(vcrit);
        ny=max(size(freqyear));
        nten=10;
        neset=max(size(eset));
        nmodel=max(size(model));
        %
        % Write to netcdf
        %
        fileout=char(strcat(file,'.nc'));
        fileoutzip=char(strcat(file,'.nc.zip'));
        %
        cmode = netcdf_getConstant('CLOBBER');
        ncid=netcdf_create(char(fileout),cmode);
        dimenn = netcdf_defDim(ncid,'Storm number',n);
        dimenm = netcdf_defDim(ncid,'Point number within storm',m);
        dimeny = netcdf_defDim(ncid,'Number of years in data set',ny);
        dimeneset = netcdf_defDim(ncid,'Eset',neset);
        dimenmodel = netcdf_defDim(ncid,'Model',nmodel);
        dimenbasin=netcdf_defDim(ncid,'Basin',2);
        dimengmeth=netcdf_defDim(ncid,'Gmeth',4);
        dimenone=netcdf_defDim(ncid,'One',1);
        dimenten = netcdf_defDim(ncid,'Ten',nten);
        %
        %  Add control variable definitions, units
        %
        basin_ID=netcdf_defVar(ncid,'Bas','NC_CHAR',dimenbasin);
        netcdf_putAtt(ncid,basin_ID,'Description','Genesis Basin: AL, AP, EP, WP, IO, SH, or GB')
        eset_ID=netcdf_defVar(ncid,'Eset','NC_CHAR',dimeneset);
        netcdf_putAtt(ncid,eset_ID,'Description','Event set name')
        fyears_ID=netcdf_defVar(ncid,'Freqyear','NC_DOUBLE',dimeny);
        netcdf_putAtt(ncid,fyears_ID,'Description','Annual frequency of events')
        freq_ID=netcdf_defVar(ncid,'Freq','NC_DOUBLE',dimenone);
        netcdf_putAtt(ncid,freq_ID,'Description','Total event frequency')
        gmeth_ID=netcdf_defVar(ncid,'Gmeth','NC_CHAR',dimengmeth);
        netcdf_putAtt(ncid,gmeth_ID,'Description','Genesis Method: rand=random seeding; clim=climatological')
        model_ID=netcdf_defVar(ncid,'Model','NC_CHAR',dimenmodel);
        netcdf_putAtt(ncid,model_ID,'Description','Name of model or reanalysis')
        vcrit_ID=netcdf_defVar(ncid,'Vcrit','NC_DOUBLE',dimenone);
        netcdf_putAtt(ncid,vcrit_ID,'Description','Minimum wind speed to be included in set')
        netcdf_putAtt(ncid,vcrit_ID,'Units','Knots')
        %
        if strcmpi(shape,'circ') && city_radius < 2001
            clat_ID=netcdf_defVar(ncid,'Clat','NC_DOUBLE',dimenone);
            netcdf_putAtt(ncid,clat_ID,'Description','Latitude of circle center')
            netcdf_putAtt(ncid,clat_ID,'Units','Degrees')
            clong_ID=netcdf_defVar(ncid,'Clong','NC_DOUBLE',dimenone);
            netcdf_putAtt(ncid,clong_ID,'Description','Longitude of circle center')
            netcdf_putAtt(ncid,clong_ID,'Units','Degrees')          
            crad_ID=netcdf_defVar(ncid,'City_radius','NC_DOUBLE',dimenone);
            netcdf_putAtt(ncid,crad_ID,'Description','Radius of circle')
            netcdf_putAtt(ncid,crad_ID,'Units','Km')
        elseif strcmpi(shape,'poly')
            polyfile=char(polyfile);
            polysize=max(size(polyfile));
            dimenpolyname=netcdf_defDim(ncid,'Polyfile',polysize);
            polyname_ID=netcdf_defVar(ncid,'Polyfile','NC_CHAR',dimenpolyname);
            netcdf_putAtt(ncid,polyname_ID,'Description','Name of poly filter file used')
            load poly.in
            [npoly,~]=size(poly);
            dimenpoly=netcdf_defDim(ncid,'Poly',npoly);
            x1_ID=netcdf_defVar(ncid,'x1','NC_DOUBLE',dimenpoly);
            x2_ID=netcdf_defVar(ncid,'x2','NC_DOUBLE',dimenpoly);
            y1_ID=netcdf_defVar(ncid,'y1','NC_DOUBLE',dimenpoly);
            y2_ID=netcdf_defVar(ncid,'y2','NC_DOUBLE',dimenpoly);
            dum_ID=netcdf_defVar(ncid,'dum','NC_DOUBLE',dimenpoly);
            netcdf_putAtt(ncid,x1_ID,'Description','Longitudes of start points of line segments')
            netcdf_putAtt(ncid,x2_ID,'Description','Longitudes of end points of line segments')
            netcdf_putAtt(ncid,y1_ID,'Description','Latitudes of start points of line segments')
            netcdf_putAtt(ncid,y2_ID,'Description','Latitudes of end points of line segments')  
            netcdf_putAtt(ncid,dum_ID,'Description','Dummy variable equal to y2')
            netcdf_putAtt(ncid,x1_ID,'Units','Degrees')
            netcdf_putAtt(ncid,x2_ID,'Units','Degrees')
            netcdf_putAtt(ncid,y1_ID,'Units','Degrees')
            netcdf_putAtt(ncid,y2_ID,'Units','Degrees')            
        end    
        %
        year_ID=netcdf_defVar(ncid,'Year','NC_DOUBLE',dimenn);
        month_ID=netcdf_defVar(ncid,'Month','NC_DOUBLE',[dimenn,dimenm]);
        day_ID=netcdf_defVar(ncid,'Day','NC_DOUBLE',[dimenn,dimenm]);
        hour_ID=netcdf_defVar(ncid,'Hour','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,hour_ID,'Units','GMT')
        lat_ID=netcdf_defVar(ncid,'Latitude','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,lat_ID,'Units','Degrees')
        lon_ID=netcdf_defVar(ncid,'Longitude','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,lon_ID,'Units','Degrees east of Greenwich')
        vc_ID=netcdf_defVar(ncid,'Vc','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,vc_ID,'Units','Knots')
        netcdf_putAtt(ncid,vc_ID,'Description','Maximum azimuthally averaged azimuthal wind')
        vnet_ID=netcdf_defVar(ncid,'Vnet','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,vnet_ID,'Units','Knots')
        netcdf_putAtt(ncid,vnet_ID,'Description', ... 
            'Maximum 10m 1 minute azimuthal wind, including effects of translation and shear')
        vinc_ID=netcdf_defVar(ncid,'Vinc','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,vinc_ID,'Units','Knots')
        netcdf_putAtt(ncid,vinc_ID,'Description','Increment to storm meridional wind to account for translation and shear')
        uinc_ID=netcdf_defVar(ncid,'Uinc','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,uinc_ID,'Units','Knots')
        netcdf_putAtt(ncid,uinc_ID,'Description','Increment to storm zonal wind to account for translation and shear')
        rm_ID=netcdf_defVar(ncid,'Rm','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,rm_ID,'Units','Km')
        netcdf_putAtt(ncid,rm_ID,'Description','Radius of maximum wind')
        router_ID=netcdf_defVar(ncid,'Router','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,router_ID,'Units','Km')
        netcdf_putAtt(ncid,router_ID,'Description','Outer radius')
        vcs_ID=netcdf_defVar(ncid,'Vcse','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,vcs_ID,'Units','Knots')
        netcdf_putAtt(ncid,vcs_ID,'Description','Maximum azimuthally averaged azimuthal wind in secondary eyewall')
        vnets_ID=netcdf_defVar(ncid,'Vnetse','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,vnets_ID,'Units','Knots')
        netcdf_putAtt(ncid,vnets_ID,'Description', ... 
            'Maximum 10m 1 minute azimuthal wind in secondary eyewall, including effects of translation and shear')
        rms_ID=netcdf_defVar(ncid,'Rmse','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,rms_ID,'Units','Km')
        netcdf_putAtt(ncid,rms_ID,'Description','Radius of secondary eyewall')
        pc_ID=netcdf_defVar(ncid,'Pc','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,pc_ID,'Units','hPa')
        netcdf_putAtt(ncid,pc_ID,'Description','Central surface pressure')
        shear_ID=netcdf_defVar(ncid,'Shear','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,shear_ID,'Units','Meters per second') 
        netcdf_putAtt(ncid,shear_ID,'Description','Environmental 850-250 hPa wind shear at storm center')
        vp_ID=netcdf_defVar(ncid,'Vp','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,vp_ID,'Units','Knots') 
        netcdf_putAtt(ncid,vp_ID,'Description','Potential intensity')
        u850_ID=netcdf_defVar(ncid,'U850','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,u850_ID,'Units','Knots') 
        netcdf_putAtt(ncid,u850_ID,'Description','Environmental zonal 850 hPa wind at storm center')
        v850_ID=netcdf_defVar(ncid,'V850','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,v850_ID,'Units','Knots') 
        netcdf_putAtt(ncid,v850_ID,'Description','Environmental meridional 850 hPa wind at storm center')
        T600_ID=netcdf_defVar(ncid,'T600','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,T600_ID,'Units','K') 
        netcdf_putAtt(ncid,T600_ID,'Description','Environmental 600 hPa temperature at storm center')
        rh600_ID=netcdf_defVar(ncid,'RH600','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,rh600_ID,'Units','percent') 
        netcdf_putAtt(ncid,rh600_ID,'Description','Environmental 600 hPa relative humidity at storm center')
        ut_ID=netcdf_defVar(ncid,'Ut','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,ut_ID,'Units','Knots') 
        netcdf_putAtt(ncid,ut_ID,'Description','Translation speed in zonal direction')
        vt_ID=netcdf_defVar(ncid,'Vt','NC_DOUBLE',[dimenn,dimenm]);
        netcdf_putAtt(ncid,vt_ID,'Units','Knots') 
        netcdf_putAtt(ncid,vt_ID,'Description','Translation speed in meridional direction')
        %
        if strcmpi(shape,'poly')
            jint_ID=netcdf_defVar(ncid,'Jint','NC_INT',[dimenn,dimenten]);
            netcdf_putAtt(ncid,jint_ID,'Description','Number index of line segment crossed')
            kint_ID=netcdf_defVar(ncid,'Kint','NC_INT',[dimenn,dimenten]);
            netcdf_putAtt(ncid,kint_ID,'Description','Along track index of line segment crossings')
            kfrac_ID=netcdf_defVar(ncid,'Kfrac','NC_DOUBLE',[dimenn,dimenten]);
            netcdf_putAtt(ncid,kfrac_ID,'Description','kfrac is the fraction of the distance between track points i and i+1 that the intersection occurs')
            xint_ID=netcdf_defVar(ncid,'Xint','NC_DOUBLE',[dimenn,dimenten]);
            netcdf_putAtt(ncid,xint_ID,'Description','Longitudes of intersections')
            netcdf_putAtt(ncid,xint_ID,'Units','Degrees')
            yint_ID=netcdf_defVar(ncid,'Yint','NC_DOUBLE',[dimenn,dimenten]);
            netcdf_putAtt(ncid,yint_ID,'Description','Latitudes of intersections')
            netcdf_putAtt(ncid,yint_ID,'Units','Degrees')
        end    
        %
        %  Compress files (decreases size of .nc files but increases the size of zipped files)
        %{
        complev=5;
        netcdf_defVarDeflate(ncid,year_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,month_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,day_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,hour_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,lat_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,lon_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,vc_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,vnet_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,rm_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,vcs_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,vnets_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,rms_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,pc_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,shear_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,vp_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,u850_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,v850_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,T600_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,rh600_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,ut_ID,true,true,complev);
        netcdf_defVarDeflate(ncid,vt_ID,true,true,complev);
        %}
        netcdf_endDef(ncid);
        %
        netcdf_putVar(ncid,eset_ID,char(eset));    
        netcdf_putVar(ncid,model_ID,char(model));        
        netcdf_putVar(ncid,basin_ID,bas);
        netcdf_putVar(ncid,gmeth_ID,gmeth);
        netcdf_putVar(ncid,vcrit_ID,vcrit);
        %
        netcdf_putVar(ncid,freq_ID,freq);        
        netcdf_putVar(ncid,fyears_ID,freqyear);
        %       
        netcdf_putVar(ncid,year_ID,yearstore);
        netcdf_putVar(ncid,month_ID,monthstore);
        netcdf_putVar(ncid,day_ID,daystore);
        netcdf_putVar(ncid,hour_ID,hourstore);
        netcdf_putVar(ncid,lat_ID,latstore);
        netcdf_putVar(ncid,lon_ID,longstore);
        netcdf_putVar(ncid,vc_ID,vstore);
        netcdf_putVar(ncid,vnet_ID,vnet);
        netcdf_putVar(ncid,rm_ID,rmstore);
        netcdf_putVar(ncid,router_ID,router);
        netcdf_putVar(ncid,vcs_ID,vsestore);
        netcdf_putVar(ncid,vnets_ID,vnetse);
        netcdf_putVar(ncid,rms_ID,rmsestore);
        netcdf_putVar(ncid,vinc_ID,vinc);
        netcdf_putVar(ncid,uinc_ID,uinc);
        netcdf_putVar(ncid,pc_ID,pstore);
        netcdf_putVar(ncid,shear_ID,shearstore);
        netcdf_putVar(ncid,vp_ID,vpstore);
        netcdf_putVar(ncid,u850_ID,u850store);
        netcdf_putVar(ncid,v850_ID,v850store);
        netcdf_putVar(ncid,T600_ID,T600store);
        netcdf_putVar(ncid,rh600_ID,rhstore);   
        netcdf_putVar(ncid,ut_ID,ut)
        netcdf_putVar(ncid,vt_ID,vt)
        if strcmpi(shape,'circ') && city_radius < 2001
            netcdf_putVar(ncid,crad_ID,city_radius)
            netcdf_putVar(ncid,clat_ID,clat)
            netcdf_putVar(ncid,clong_ID,clong)
         elseif strcmpi(shape,'poly')
             netcdf_putVar(ncid,polyname_ID,polyfile)
             netcdf_putVar(ncid,jint_ID,jint)
             netcdf_putVar(ncid,kint_ID,kint)
             netcdf_putVar(ncid,kfrac_ID,kfrac)
             netcdf_putVar(ncid,xint_ID,xint)
             netcdf_putVar(ncid,yint_ID,yint)
             load poly.in
             x1=poly(:,1);y1=poly(:,2);x2=poly(:,3);y2=poly(:,4);
             %disp(max(y2))
             netcdf_putVar(ncid,x1_ID,double(x1))
             netcdf_putVar(ncid,x2_ID,double(x2))
             netcdf_putVar(ncid,y1_ID,double(y1))
             netcdf_putVar(ncid,y2_ID,double(y2))       
             dum=y2;
             netcdf_putVar(ncid,dum_ID,double(dum))
        end 
        %
        netcdf_close(ncid) 
        %zip(fileoutzip,fileout)
        %pause(3)
        %       
