close all


% Read SSJ4 .txt file for a day 
hemisphere = 'N';
hemisphere = 'S';
%filname = '/home/liamk/data/cenergy/cenergy_all_f15_20050514.txt';
%filname = '/home/emery/reu11/cenergy_all_f15_20050110.txt';
nskp = 70801;
nrd = 72587-nskp;
poymax = 20.028;
poyx = [20.0 20.13];
emax = 17.;
%filname = '/home/emery/reu11/cenergy_all_f15_20041107.txt';
%poyname = '/hao/aim3/dmsp/poynting_flux/pfomni_04312.txt';
% 19UT 67763-69404 (fracday mlat MLT ? ? ? PoynF ? ? By Bz Vsw Dsw ? of 14)
nskp = 67211;
nrd = 68965-nskp;
nskpp = 67763
nrdp = 69404-nskpp;
%poymax = 19.15;
%poyx = [19.1 19.3];
emax = 10.;
poyn = -200.;
%filname = '/home/emery/reu11/cenergy_all_f15_20001108.txt';
filname = '/hao/aim3/dmsp/dmsp_rich/txt_dmspssj/2005/13_20050116.txt';
% 22104-23648 for SH 6.3UT minpoy=-0.03!!  (21549-23226ssj)  and nn=0 problems!
hemisphere = 'S';
nskp = 21549;
nrd = 14851-nskp;
nskpp = 22104;
nrdp = 23648-nskpp;
x1 = 6.0;
x2 = 7.0;
x3 = 0;
x4 = 24;
poyn = -50.;
%12914-14504 for NH 3.7UT (13156-14851ssj4) 
hemisphere = 'S';
nskp = 13156;
nrd = 14851-nskp;
nskpp = 12914;
nrdp = 14504-nskpp;
x1 = 3.6;
x2 = 4.2;
emax = 5;
poyn = -12.;
kmax = 12.;
fxmin = 0.03;
fxmin = 0.04;
cos_85 = 0.087155743;
Sa = 128.30;
badmlat = 70.;
detrendmlat = 60.; 
mmax = 24; 
nsm = 25;
bdpt = 0;
eflxmin = 0.3;
tossx = 20.;

%Getting the date from the filename

inst   = filname(43:44);
iyear  = filname(46:49);
iyr    = filname(48:49);  
iday   = filname(52:53);
imonth = filname(50:51);

PlotTitle = ([imonth, ' ', iyear,' ',iday,' ',hemisphere,'H DMSP-F', inst,'SSJ4'] );

iyearN = str2num(iyear);
idayN  = str2num(iday);

% Read Poynting day 1 sec file with OMNI data

% Read SSJ4 data
SSJ=fopen(filname,'r');

for n=1:nskp
 topssj = fgetl(SSJ);
end
nn = 0;
for n=1:nrd
 lindat = fgetl(SSJ);
 weidat = sscanf(lindat,'%f');
 if abs(weidat(8)) > 45
  nn = nn + 1;
  uts(nn)  = weidat(2);
  glat(nn) = weidat(6); 
  glon(nn) = weidat(7);
  mlat(nn) = weidat(8);
  mlon(nn) = weidat(9);
  alat(nn) = abs(weidat(8));
  mlt(nn)  = weidat(10);
  eflx(nn) = weidat(11);
  ekev(nn) = weidat(13)*0.001;
  iflx(nn) = weidat(14);
  ikev(nn) = weidat(16)*0.001;
  if (eflx(nn)<fxmin)
   ekev(nn) = 0.;
  end
  if (iflx(nn)<fxmin)
   ikev(nn) = 0.;
  end
 end
end
uts1 = uts(1);
% 05134 nn=1675 but uts is filled to 21264 values for some reason.  Therefore limit plot to nn
utslen  = length(uts);
mlatlen = length(mlat);
minuts  = min(uts(1:nn));
maxuts  = max(uts(1:nn));
maxeflx = max(eflx(1:nn));
maxiflx = max(iflx(1:nn));
maxekev = max(ekev(1:nn));
maxikev = max(ikev(1:nn));
minmlat = min(mlat(1:nn));
emax    = max(eflx)+0.5;


% Conductivity Calculations
% Calculating Particle Induced Conductivity 

for n = 1:nn
	condp(n) = 40 * (ekev(n)) * (eflx(n))^(1/2) / (16 + (ekev(n))^2);

	condh(n) = 0.45 * (ekev(n))^0.85 * condp(n);

end 

% Calculating conductivities from UV/EUV
% Calls in variables from subsol.m and cossza.m
%   csza = solar zenith angle
%   slat = co-latitude of the sub solar point
%   slon = longitude of the sub solar point
%   Sa   = solar radio flux in W/m^2 Hz
%   condUVp = conductivity from solar UV, EUV, and cosmic rays
%   condp_Tot = Total Pederson conductivity
	
 for n = 1:nn
	
	[slat(n) slong(n)] = subsol(iyearN, idayN, uts(n)/3600, uts(n)/60, uts(n));

	csza(n) = cossza(glat(n), glon(n), slat(n), slong(n));	

	if csza(n) > cos_85

  		condUVp(n) = 0.88 * (Sa * csza(n))^(1/2) + 0.1; 
		condUVh(n) = 1.5 * (Sa * csza(n))^(1/2) +0.2;

	else 
		f(n)  = 0.06 + exp(1.803 * tanh(3.833 * csza(n)) + (0.5 * csza(n)) - 2.332);  
		
		Bp(n) = (1 - (0.99524 * (sin(slat(n)))^2)) * (1 + (0.3 * cos(slat(n))^2)); 

		condUVp(n) = 12.5 * (Sa/180)^(1/2) * f(n) / Bp(n);
	end	
	
	condp_Tot(n) = (condp(n)^2 + condUVp(n)^2)^(1/2);		
end	


%----------------------Tosses sharp peaks in Eflux data-------------------------
% Tossing peaks greater than 20x unsmoothed mean
% Prints out number of points tossed

meanflx = mean(eflx);

for n = 1:utslen

	if eflx(n) < tossx*meanflx
		mflx(n) = eflx(n);
	else 
		mflx(n) = NaN;
		bdpt = bdpt+1;
	end
end 

bdpt

% Linearly interpolates gaps in time series 
% R. Pawlowicz 6/Nov/99

mflx2=mflx;

bd=isnan(mflx);
gd=find(~bd);

bd([1:(min(gd)-1) (max(gd)+1):end])=0;


mflx2(bd)=interp1(gd,mflx(gd),find(bd));

% Smooth(y) does a moving smoothing over default of 5 values (20 sec for 4 sec data); 
% Smooth(y,span) where span MUST BE ODD - hence smooth(y,10) is converted to the default 5
% Smooths the new eflx set of values with the tossed peaks > 8x mean 
% Smooths the mean energy <E> KeV 
% Smooths condh data

aveflx = smooth( mflx2,nsm); 
mekev  = smooth( ekev, nsm);
avecondh = smooth(condh,nsm);

% Filtering bad data on Ekev plot
% Do not plot E-kev if E-flux less than 0.20

for  n = 1:nn 

	if aveflx(n) <= 0.20

		ekev(n)  = 0;
		mekev(n) = 0;
	end
end 

% Will label start of pass in hours=uts/3600.
% data1 = [uts/3600. alat mlt eflx];  
for n=1:nn
data2(n,1) = uts(n)/3600.;
data2(n,2) = mlat(n);
data2(n,3) = mlt(n);
data2(n,4) = eflx(n);
end
%data2 = [uthp(:,1) mlatp(:,1) mltp(:,1) abs(poyf(:,1))];	% works!
   [out,track,lines] = satelliteTrackPlot(data2,[0 10], 15, PlotTitle, 1,0,{});

%stop

% ------------------Finding Peaks and 1/e Fall Off in Eflux------------------

% Using the index of the maximum mlat (ixmlat) to find the maximum sunward velocities on either side of this index

xmlat = 0.5*(badmlat+detrendmlat);
ixmlat = 0;
for i=1:utslen
  if (abs(mlat(i)) > xmlat)
    xmlat = abs(mlat(i));
    ixmlat = i;
  end
end

% Use the max funtion to find the two peaks

[peak1 ipeak1] = max(aveflx(1:ixmlat));

[peak2 ipeak2] = max(aveflx(ixmlat:mlatlen));

ipeak2 = ipeak2 + ixmlat; 

peakarr = [peak1 peak2];
ipeakarr = [ipeak1 ipeak2];

np = length(peakarr);

% Calculate 1/e fall off for each peak

peakfall = peakarr(1:np)/2.718;

% Finding the 1/e fall off value before and after each peak (4 total values) closest to equator and poles

for i = 1:ipeak1

	if aveflx(i) > peakfall(1) 
		ifall1 = i-1;  
		break
	end
end


for i = ixmlat:-1:ipeak1

	if aveflx(i) > peakfall(1) 
		ifall2 = i-1;  
		break
	end
end


for i = ixmlat:ipeak2

	if aveflx(i) > peakfall(2) 
		ifall3 = i-1;  
		break
	end
end

for i = mlatlen:-1:ipeak2

	if aveflx(i) > peakfall(2) 
		ifall4 = i;  
		break
	end
end

fall = [ifall1 ifall2 ifall3 ifall4];

% Finding endpoints for Gaussian curve (appx. where peak approaches zero)

for i = 1:utslen
	if i < fall(1)	
		if aveflx(i) < peakarr(1)*.05 && aveflx(i) 
		istart(1) = i;
		start(1) = aveflx(i);
		end
	elseif i >= fall(2) && i < ixmlat
		[start(2) istart(2)] = min(aveflx(ifall2:ixmlat));
		istart(2) = istart(2) + fall(2);
	elseif i >= ixmlat && i < fall(3)
		[start(3) istart(3)] = min(aveflx(ixmlat:ifall3));
		istart(3) = istart(3) + ixmlat;
	elseif i >= fall(4) && i <= utslen
		if aveflx(i) < peakarr(2)*.05 
		istart(4) = i;
		start(4) = aveflx(i)
		break
		end
	end
end		


% Calculates the Gaussian curves 
% Calculates the summation of each of the four half-widths from both the Gaussian and original EFlux curve
% Calculated the MLAT delta h

gaussum = [0 0 0 0];
origsum = [0 0 0 0];

for i = istart(1):ipeak1
	
	gaussflx1(i) = peakarr(1) * exp( -(abs(mlat(ipeakarr(1))-mlat(i)))^2 / (mlat(ipeakarr(1))-mlat(fall(1)))^2 );
	gaussum(1)   = gaussum(1) + gaussflx1(i);
	origsum(1)   = origsum(1) + aveflx(i);
	delh(1)      = abs(mlat(ipeakarr(1))-mlat(fall(1)));
end

for i = ipeak1:istart(2)

	gaussflx2(i) = peakarr(1) * exp( -(abs(mlat(ipeakarr(1))-mlat(i)))^2 / abs(mlat(ipeakarr(1))-mlat(fall(2)))^2 );
	gaussum(2)   = gaussum(2) + gaussflx2(i);
	origsum(2)   = origsum(2) + aveflx(i);
	delh(2)      = abs(mlat(ipeakarr(1))-mlat(fall(2)));
end

for i = istart(3):ipeak2

	gaussflx3(i) = peakarr(2) * exp( -(abs(mlat(ipeakarr(2))-mlat(i)))^2 / abs(mlat(ipeakarr(2))-mlat(fall(3)))^2 );
	gaussum(3)   = gaussum(3) + gaussflx3(i);
	origsum(3)   = origsum(3) + aveflx(i);
	delh(3)      = abs(mlat(ipeakarr(2))-mlat(fall(3)));
end

for i = ipeak2:istart(4)

	gaussflx4(i) = peakarr(2) * exp( -(abs(mlat(ipeakarr(2))-mlat(i)))^2 / abs(mlat(ipeakarr(2))-mlat(fall(4)))^2 );
	gaussum(4)   = gaussum(4) + gaussflx4(i);
	origsum(4)   = origsum(4) + aveflx(i);
	delh(4)      = abs(mlat(ipeakarr(2))-mlat(fall(4)));
end

delmax=max(delh(1:4))+.5;

% Calculate Original sum to Gaussian sum ratio

ratiosum(1:4) = gaussum ./ origsum;

% Getting poleward/equatorward and 06MLT/18MLT values for uts, mlat, mlt, <E>, eflx


for i = ifall1:ipeak1
	if mlt(i) >= 0 && mlt(i) < 12
		utse6(i)    = uts(i);
		avflxe6(i)  = aveflx(i); 
		mlte6(i)    = mlt(i);
		mlate6(i)   = mlat(i);
		ekeve6(i)   = mekev(i);	
	elseif mlt(i) >= 12 && mlt(i) <= 24
		uts12(i)    = uts(i);
		avflxe12(i) = aveflx(i);
		mlte12(i)   = mlt(i);
		mlate12(i)  = mlat(i);
		ekeve12(i)  = mekev(i);
	end
end

for i = ipeak1:ifall2

	if mlt(i) >=0 && mlt(i) < 12
		utsp6(i)   = uts(i);
		avflxp6(i) = aveflx(i);
		mltp6(i)   = mlt(i);
		mlatp6(i)  = mlat(i);
		ekevp6(i)  = mekev(i);
	elseif mlt(i) >= 12 && mlt(i) <= 24
		utsp12(i)  = uts(i);
		avflxp12(i)= aveflx(i);
		mltp12(i)  = mlt(i);
		mlatp12(i) = mlat(i);
		ekevp12(i) = mekev(i);
	end
end

for i = ifall3:ipeak2	

	if mlt(i) >=0 && mlt(i) < 12
		utsp6(i)   = uts(i);
		avflxp6(i) = aveflx(i);
		mltp6(i)   = mlt(i);
		mlatp6(i)  = mlat(i);
		ekevp6(i)  = mekev(i);
	elseif mlt(i) >= 12 && mlt(i) <= 24
		utsp12(i)  = uts(i);
		avflxp12(i)= aveflx(i);
		mltp12(i)  = mlt(i);
		mlatp12(i) = mlat(i);
		ekevp12(i) = mekev(i);
	end
end

for i = ipeak2:ifall4

	if mlt(i) >=0 && mlt(i) < 12
		utse6(i)   = uts(i);
		avflxe6(i) = aveflx(i);
		mlte6(i)   = mlt(i);
		mlate6(i)  = mlat(i);
		ekeve6(i)  = mekev(i);
	elseif mlt(i) >= 12 && mlt(i) <= 24
		utse12(i)  = uts(i);
		avflxe12(i)= aveflx(i);
		mlte12(i)  = mlt(i);
		mlate12(i) = mlat(i);
		ekeve12(i) = mekev(i);
	end
end


%------------------------------ Plots -----------------------------------

% top plot
% Plots E flux smoothed over 19 seconds (red), E flux unsmoothed raw data (cyan), and Gaussian curves (black) if peak>0.3

dark_green = [0 0.5 0];
fsize = 12;
fsiz = 11;
	H1=figure;
		subplot ('Position',[0.1,0.70,0.83,0.20])
	hold on;
	plot(uts(1:nn)/3600.,eflx(1:nn),'c');
	plot(uts(1:nn)/3600.,aveflx(1:nn),'m');
	plot(uts(istart(1:4))/3600.,start(1:4),'k.');
if peakarr(1) > eflxmin
	plot(uts(istart(1):ipeak1)/3600.,gaussflx1(istart(1):ipeak1),'k', 'LineWidth',2);
	plot(uts(ipeak1:istart(2))/3600.,gaussflx2(ipeak1:istart(2)),'k', 'LineWidth',2);
end

if peakarr(2) > eflxmin
	plot(uts(istart(3):ipeak2)/3600.,gaussflx3(istart(3):ipeak2),'k', 'LineWidth',2);
	plot(uts(ipeak2:istart(4))/3600.,gaussflx4(ipeak2:istart(4)),'k', 'LineWidth',2);
end
	text (uts(1)/3600.,0.65*emax,'e-','Color','c','fontsize',fsize, 'fontweight','bold');
	text (uts(1)/3600.,0.45*emax,'. = 1/e','color','m','fontsize',fsize, 'fontweight','bold');
	text (uts(1)/3600.,0.25*emax,'\diamondsuit = Peak','color','r','fontsize',fsize, 'fontweight','bold');

% Plotting peaks and 1/e fall offs on top of eflx plot

for i = 1:np
	
	if peakarr(i) > eflxmin
		plot(uts(ipeakarr(i))/3600, peakarr(i), 'rd', 'MarkerFaceColor', 'r');	
	end 
end
	
if peakarr(1) > eflxmin 	
	plot(uts(fall(1:2))/3600, peakfall(1), 'm.');
end
	
if peakarr(2) > eflxmin
	plot(uts(fall(3:4))/3600, peakfall(2), 'm.');	
end	

grid;
                        ylabel('E flux (mW/m^2)','fontsize',fsize,'fontweight','bold')
	                title (PlotTitle,'fontsize',fsize,'fontweight','bold');
	set (gca,'Xlim',[x1 x2],'Ylim',[0 emax],'XTickLabel',{' ';});
		subplot ('Position',[0.1,0.50,0.83,0.20])
hold on;
	plot(uts(1:nn)/3600.,ekev(1:nn),'g');
	plot(uts(1:nn)/3600.,mekev(1:nn),'k');

% MLT and Dsw plots to save space
%	plot(uts(1:nn)/3600.,mlt(1:nn),'b');
%       plot(uthp,dsw,'k');
grid;
         ylabel('<E> (keV)','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[0 kmax],'XTickLabel',{' ';});

    
% Plotting conductivities
subplot ('Position',[0.1,0.30,0.83,0.10])
hold on;
	
	plot(uts(1:nn)/3600., condp(1:nn),'r'); 
%	plot(uts(1:nn)/3600., avecondh(1:nn), 'y');
	plot(uts(1:nn)/3600., condp_Tot(1:nn), 'k');
	plot(uts(1:nn)/3600., condUVp(1:nn), 'm');
	text (uts(1)/3600.,10,'\Sigma_P','Color','r','fontsize',fsize, 'fontweight','bold');
	text (uts(100)/3600.,10,'\Sigma_PTot','color','k','fontsize',fsize, 'fontweight','bold');
	text (uts(300)/3600.,10,'\Sigma_PUV/EUV','color','m','fontsize',fsize, 'fontweight','bold');
        
grid;
        ylabel('\Sigma_P, \Sigma_PTot', 'fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[0 20],'XTickLabel',{' ';});

% bottom plots MLAT and MLT
% Plot MLT and Dsw 
		subplot ('Position',[0.1,0.20,0.83,0.10])
	hold on;
	plot(uts(1:nn)/3600.,mlt(1:nn),'b');

grid;
        ylabel('MLT','fontsize',fsiz,'fontweight','bold')
	set (gca,'Xlim',[x1 x2],'Ylim',[0 24],'XTickLabel',{' ';});
% Plot MLAT axis
		subplot ('Position',[0.1,0.10,0.83,0.10])
	hold on;
	plot(uts(1:nn)/3600.,mlat(1:nn),'b');

grid;
        ylabel('MLAT','fontsize',fsiz,'fontweight','bold')
        xlabel('UT (h)','fontsize',fsize,'fontweight','bold')
	set (gca,'Xlim',[x1 x2]);


hold off;

% Figure (3) 
% First(Top) Plot: plots the ratio of the E Flux fitted Gaussian curve over the raw data curve
% Second Plot: plots the peaks of the E Flux
% Third Plot: plots the MLAT delta h for each 1/e fall off

figure(3)

	subplot ('Position',[0.1,0.60,0.83,0.30])
hold on;
	plot(mlt(fall(2:3)),ratiosum(2:3),'c*',mlt(fall(1)),ratiosum(1),'r*', mlt(fall(4)),ratiosum(4),'r*');	
	text (2,1.8,'* Poleward fall off','Color','c','fontsize',fsize, 'fontweight','bold');
	text (2,1.5,'* Equatorward fall off','Color','r','fontsize',fsize, 'fontweight','bold');
grid;
hold on;
	plot((0:mmax), ones(1,25),'color','k');
	title ('Peak and Fall Off Comparisons','fontsize',fsize,'fontweight','bold');
	ylabel('Ratio Fit/Data','fontsize',fsiz,'fontweight','bold');
	set (gca,'Xlim',[x3 x4], 'Ylim',[0 2]);	
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	subplot ('Position',[0.1,0.40,0.83,0.20])
grid;
hold on;
	set (gca,'Xlim',[x3 x4], 'Ylim',[0 emax]);
	plot(mlt(ipeakarr(1:2)), peakarr(1:2),'md','MarkerFaceColor','m');
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');	
	text (2,.8*emax,'Peaks','Color','m','fontsize',fsize, 'fontweight','bold');
	ylabel('E Flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
	subplot ('Position',[0.1,0.10,0.83,0.30])
hold on;
	set (gca,'Xlim',[x3 x4], 'Ylim',[0,delmax]);	
	plot(mlt(fall(2:3)),delh(2:3),'c^','MarkerFaceColor','c')
hold on;
	plot(mlt(fall(1)),delh(1),'r^','MarkerFaceColor','r');
	plot(mlt(fall(4)),delh(4),'r^','MarkerFaceColor','r');
	yL = get(gca,'YLim');
	line([6 6],yL,'color','k','LineStyle','--');
	line([18 18],yL,'color','k','LineStyle','--');
	text (6.5,0.2,'6 MLT','fontsize',fsize, 'fontweight','bold');
	text (18.5,0.2,'18 MLT','fontsize',fsize, 'fontweight','bold');
	text (2,0.8*delmax,'Poleward \Deltah','Color','c','fontsize',fsize, 'fontweight','bold');
	text (2,0.7*delmax,'Equatorward \Deltah','Color','r','fontsize',fsize, 'fontweight','bold');
	ylabel('MLAT \Deltah','fontsize',fsiz,'fontweight','bold');
	xlabel('MLT','fontsize',fsiz,'fontweight','bold');	
grid;
hold off;

figure(4)
	plot(ekeve6, avflxe6, 'r+',ekevp6,avflxp6,'b+');
hold on;
	plot(ekeve12, avflxe12, 'r+',ekevp12,avflxp12,'b+');
	text (0.5,0.8*max(aveflx),'Poleward','Color','b','fontsize',fsize, 'fontweight','bold');
	text (0.5,0.7*max(aveflx),'Equatorward','Color','r','fontsize',fsize, 'fontweight','bold');
	title('E-Flux and Mean Energy Relation','fontsize',fsize,'fontweight','bold');
	ylabel('E flux (mW/m^2)','fontsize',fsiz,'fontweight','bold');
	xlabel('<E> (keV)','fontsize',fsiz,'fontweight','bold');
hold off;


% Write out data to text file:
% Year, Day, Time averages for smoothed data (nsm), Minimum e-flux (eflxmin), Number times the mean over which peaks are  tossed (tossx), UTS, Smoothed E Flux, MLAT, MLT, Mean energy <E> keV
% Values are only for indices within 1/e fall off to peak range


fout1 = sprintf('%s%s%seq06.txt',iyr,imonth,iday);
HPOUT1 = fopen( fout1, 'w' );
fprintf( HPOUT1, '%4.0d %2.0d %2.0d %.2d %.2f \n\n %5.0f %6.2f %6.2f %6.2f %6.2f \n', iyearN, idayN, nsm, eflxmin, tossx, utse6, avflxe6, mlte6, mlate6, ekeve6);

fout2 = sprintf('%s%s%seq12.txt', iyr, imonth, iday);
HPOUT2 = fopen(fout2, 'w' );
fprintf( HPOUT1, '%4.0d %2.0d %2.0d %.2d %.2f \n\n %5.0f %6.2f %6.2f %6.2f %6.2f \n', iyearN, idayN, nsm, eflxmin, tossx, utse12, avflxe12, mlte12, mlate12, ekeve12);

fout3 = sprintf('%s%s%spo06.txt', iyr, imonth, iday);
HPOUT3 = fopen(fout3, 'w');
fprintf( HPOUT1, '%4.0d %2.0d %2.0d %.2d %.2f \n\n %5.0f %6.2f %6.2f %6.2f %6.2f \n', iyearN, idayN, nsm, eflxmin, tossx, utsp6, avflxp6, mltp6, mlatp6, ekevp6);

fout4 = sprintf('%s%s%spo12.txt', iyr, imonth, iday);
HPOUT4 = fopen(fout4, 'w');
fprintf( HPOUT1, '%4.0d %2.0d %2.0d %.2d %.2f \n\n %5.0f %6.2f %6.2f %6.2f %6.2f \n', iyearN, idayN, nsm, eflxmin, tossx, utsp12, avflxp12, mltp12, mlatp12, ekevp12);




