disp('AMMA_da_red1')
%toga coare bulk flux model version 3.0
%***************************************
%uses following matlab subroutines:
%	julday.m
%	cor30a.m
%	psiu_30.m
%	psit_30.m
%	grav.m
%***************************************
%*********** basic specifications  *****
%	zu=			height of wind measurement
%	zt=			height of air temperature measurement
%	zq=			height of air humidity measurement
%	ts_depth		depth of water temperature measurement
%	jwarm=			0=no warm layer calc, 1 =do warm layer calc
%	jcool=			0=no cool skin calc, 1=do cool skin calc
%***********   input data **************
%	jdy=		date in decimal julian day format
%	lon=			longitude (deg e=+)
%	u=			wind speed (m/s), height zu
%	us=			surface current (m/s)
%	ts=			bulk surface sea temp (cent)
%	t=			air temp (cent), height zt
%	qs=			sea surface sat specific humidity (g/kg)
%	q=			air specific humidity (g/kg), height zq
%	Rs=			downward solar flux (w/m^2)
%	Rl=			downward IR flux (w/m^2)
%	zi=			inversion height (m)
%	P=			air pressure (mb)
%	rain=			rain rate (mm/hr)
%	lat=		lat
%********** output data  ***************
%	hsb=			sensible heat flux (w/m^2)
%	hlb=			latent heat flux (w/m^2) inlcuding webb correction
%	RF=			rain heat flux(w/m^2)
%	W_webb=		mean webb velocity (m/s)
%	tau=			stress (nt/m^2)
%	zo=			velocity roughness length (m)
%	zot			temperature roughness length (m)
%	zoq=			moisture roughness length (m)
%	L=			Monin_Obukhov stability length
%	usr=			turbulent friction velocity (m/s), including gustiness
%	tsr			temperature scaling parameter (K)
%	qsr			humidity scaling parameter (g/g)
%	dter=			cool skin temperature depression (K)
%	dqer=			cool skin humidity depression (g/g)
%	tkt=			cool skin thickness (m)
%	Cd=			velocity drag coefficient at zu, referenced to u
%	Ch=			heat transfer coefficient at zt
%	Ce=			moisture transfer coefficient at zq
%	Cdn_10=			10-m velocity drag coeeficient, including gustiness
%	Chn_10=			10-m heat transfer coeeficient, including gustiness
%	Cen_10=			10-m humidity transfer coeeficient, including gustiness
%
fclose('all');clear all;
%
didread=0;
if didread~=1;
    x=load('C:\Data\AMMA\flux_data\procturb\AMMA08_da.txt');%read file with half-hour average data; set your local directory 
      yxx=load('C:\Data\AMMA\flux_data\procturb\AMMA08_cr.txt');%read file with licor data
      yco=load('C:\Data\AMMA\flux_data\procturb\AMMA08_co.txt');
      yfl=load('C:\Data\AMMA\flux_data\procturb\AMMA08_fl.txt');
         
   jdy=x(:,1);%julian day at beginning of half hour average
   %*****   ships imet and scs data
   U_scs=x(:,2); %true wind speed, m/s; imet sonic anemometer
   dir_scs=x(:,34);%true wind direction (from), deg (clockwise rel north)
   tsg=x(:,4);%tsg water temperature, C (5 m depth)
   ta_im=x(:,5);%imet air temperature, C 
   qs_tsg=x(:,6);%imet sea surface specific humidity, g/kg 
   qa_im=x(:,7);%imet air specific humidity, g/kg
   rs_im=x(:,8);%downward solar flux, W/m^2 (imet units)
   rl_im=x(:,9);%downward ir flux, W/m^2 (imet units)
   sig_rl_im=x(:,157);
   rain_im=x(:,10);%
   zu_im=x(:,13);%
   zt_im=x(:,14);%
   zq_im=x(:,15);%
	%**********   ETL data  **********   
   U=x(:,16);%true wind from ETL sonic
	dir=x(:,41);%
   tsnk=x(:,18);%sea snake temperature, C (0.05 m depth)
	ta=x(:,19);%air temperature, ETL   
   qse=x(:,20);%sea surface specific humidity, g/kg (from snake)
   qa=x(:,21);%air specific humidity, g/kg (z=15  m)
   rs=x(:,22);%downward solar flux, W/m^2 (ETL units)
   sig_rl=x(:,95);%STD downward solar flux, W/m^2 (ETL units)
	rl=x(:,23);%downward IR flux, W/m^2 (ETL units)
	sig_rs=x(:,96);%STD downward IR flux, W/m^2 (ETL units)
   org=x(:,24);%rainrate, mm/hr (ETL STI optical rain gauge, uncorrected)
	zu_etl=x(:,27);%
   zt_etl=x(:,28);
   zq_etl=x(:,29);
   ustr=x(:,60);
   %************  Bulk Fluxes  ************
   Hsb=x(:,74);%bulk sensible heat flux, W/m^2
   Hlb=x(:,75);%bulk latent heat flux, W/m^2 (includes Webb et al. correction?)
   Taub=x(:,76);%bulk wind stress along mean wind, N/m^2
      
   ustr0=x(:,81);
   tstr0=x(:,82);
   qstr0=x(:,83);
   Hrain=x(:,87);%rain heat flux, W/m^2 as per Gosnell et al
   
   %************   Navigation & data quality indicators
   Lat=x(:,151);%latitude, deg  (SCS pcode)
   Lon=x(:,150);%longitude, deg (SCS pcode)
   urel=x(:,38);%relative wind speed, m/s
   reldir=x(:,42);%relative wind direction (from), deg clockwise rel ship's bow(ETL sonic)
   tiltx=x(:,44)*180/pi;%flow tilt, deg
	head=x(:,112);%ship heading, deg clockwise rel north (SCS laser ring gyro)
   sig_h=x(:,113);%std dev of heading, >5 indicates ship manauver 
   ushp_im=x(:,97);%ship speed, m/s (SCS GPS)
   sig_uim=x(:,98);%std dev of ship speed, m/s (>.2 indicates manauver during average)
   sig_him=x(:,102);%std dev of ship cog)
   ushp=x(:,153);%doppler log
   sig_u=x(:,154);
   rho=x(:,52);
   %################################comparison ludo#####################################
%    figure;plot(jdy,U,'r.',jdy,U_scs,'.');title(' True wind speed (m/s), blue=ship, red=ETL');
%    figure;plot(jdy,dir,'r.',jdy,dir_scs,'.');title(' True wind direction (m/s), blue=ship, red=ETL');
%    figure;plot(jdy,tsnk,'r.',jdy,tsg,'.');title(' SST(degC), blue=ship tsg, red=ETL seasnake');
%    figure;plot(jdy,ta,'r.',jdy,ta_im,'.');title('Air temperature (degC), blue=ship, red=ETL');
%    figure;plot(jdy,qa,'r.',jdy,qa_im,'.');title('Air specific humidity (g/kg), blue=ship, red=ETL');
%    figure;plot(jdy,rs,'r.',jdy,rs_im,'.');title('downward solar flux, W/m^2, blue=ship, red=ETL');
%    figure;plot(jdy,rl,'r.',jdy,rl_im,'.');title('downward ir flux, W/m^2, blue=ship, red=ETL');
%    figure;plot(jdy,org,'r.',jdy,rain_im,'.');title('rainrate, mm/hr, blue=ship, red=ETL');
   %*************   Turbulence Data  ***********
%    wu1=x(:,116);%method 1
%    wu=x(:,117);
%    wv=x(:,119);
%    wt_son=x(:,121);
%    wq=x(:,123);
%    uu=x(:,127);
%    vv=x(:,128);
%    ww=x(:,129);
%    tt=x(:,130);
%    qq=x(:,131);
   wu=x(:,116);%method 1
   wuj=x(:,117);%edson
   wv=x(:,118);
   wvj=x(:,119);%edson
   wt_son=x(:,121);
   %wq=x(:,123);
   uu=x(:,122);
   uuj=x(:,125);
   vvj=x(:,126);
   vv=x(:,123);
   wwj=x(:,127);
   ww=x(:,124);
   tt=x(:,128);
   %qq=x(:,131)
   
   %sigoph=x(:,133);
   %ophm=x(:,134);
   
   sp1=x(:,129);%motion corrections, u (m/s)
   sp2=x(:,130);%motion corrections, v (m/s)
   sp3=x(:,131);%motion corrections, w (m/s)
   
   cu=x(:,135);%cu^2
   cw=x(:,137);%cw^2
   ct=x(:,139);%ct^2
   co=x(:,141);%cq^2 licor???
   cql=yco(:,243); %cq^2 Licor
%    figure;plot(cq,'r');hold on;plot(cql);axis([0 2500 0 0.1]);
%    title('Cq^2 form Ophir (red) and Licor (blue)');
   
   usi=x(:,143);%inertial dissipation scaling paramters
   tsi=x(:,145);
   osi=x(:,147);
  %*********  fast humidity sensors 
   %Q_oph=x(:,164); %absolute humidity from Ophir unit (g/m3)
   %wq_oph=x(:,123);
   %sgq_oph=sqrt(x(:,131));
   %ii=find(sgq_oph>1);
   %sgq_oph(ii)=NaN;Q_oph(ii)=NaN; wq_oph(ii)=NaN;
   
   q_lic=yxx(:,2); %specific humidity from Licor unit (g/kg)
   wq_lic=yxx(:,9);
   wco2_lic=yxx(:,11);
   P_lic=yxx(:,6);%internal pressure on Licor
   sgP_lic=yxx(:,16);%
   sgq_lic=yxx(:,12);
    
   co2_lic=yxx(:,4);
   ii=find(co2_lic<350 | co2_lic>400);co2_lic(ii)=NaN;
   yxx(ii,3:6)=NaN;yxx(ii,10:13)=NaN;
   
   sgq_lic=yxx(:,12);
   ii1=find(sgq_lic>1);
   yxx(ii1,3:6)=NaN;yxx(ii1,10:13)=NaN;
   sgq_lic(ii1)=NaN;
    
   sgco2_lic=yxx(:,14);
   ii2=find(sgco2_lic>4);
   yxx(ii2,3:6)=NaN;yxx(ii2,10:13)=NaN;
   sgco2_lic(ii2)=NaN;

   
q_lic=yxx(:,2);
wq_lic=yxx(:,9);
wco2_lic=yxx(:,11);
co2_lic=yxx(:,4);
lagtime=yfl(:,2);   %10min mean lagtime
O3mean_o=yfl(:,3);  %crosscorrelation technique
O3dev_o=yfl(:,4);   %stdev O3 concentration
O3w_o=yfl(:,5);     %fluxes calculated by crosscorrelation and 10 min lagtime
O3w_h=yfl(:,6);     %fluxes calculated by crosscorrelation and hourly (average 6x10 min mean) lagtime
O3w_d=yfl(:,7);     %fluxes calculated by crosscorrelation and dayly lagtime (average day, 24x6x10 minmean)
O3mean=yfl(:,8);    % first 51 values dumped
O3dev=yfl(:,9);     %stdev O3 concentration
O3w=yfl(:,10);      %fluxes calculated with constant lagtime of 5.1 s
slowO3=yfl(:,11);      %O3 conc Ml8810
slowO3std=yfl(:,12);

figure;plot(jdy,O3mean,'r.',jdy,O3mean_o,'.');title('O3 / ppb, blue=crosscorrel, red=const lagtime');xlabel('Day of Year');ylabel('O3 / ppb');
figure;plot(jdy,O3w,'r.',jdy,O3w_o,'.');title('wO3 / (cm/s.)ppb, blue=crosscorrel, red=const lagtime');xlabel('Day of Year');ylabel('wO3 / (cm/s.)ppb');
figure;plot(jdy,O3w_o,'.',jdy,O3w_h,'r.',jdy,O3w_d,'g.');title('wO3 / (cm/s.)ppb, crosscorrel, blue= tenmin lagtime, red=hourly lagtime, green=daily lagtime+limits');xlabel('Day of Year');ylabel('wO3 / (cm/s.)ppb')


end;


%*******************  set constants  ****************
	tdk=273.16;
	grav=9.72;
	Rgas=287.1;
	cpa=1004.67;
	rhoair=1010*100./(Rgas*(ta+tdk).*(1+0.61*qa/1000));
	ii=find(isfinite(tsnk));Le=(2.501-.00237*median(tsnk(ii)))*1e6;
	
	be=0.026;
	cpw=4000;
	rhow=1022;
	visw=1e-6;
	tcw=0.6;
   
%***********************   make modifications, calibrations, limit tests, etc
AMMA_fixit_da1;

%***********************   table to account for sea snake out of water
tm_in=[192 209.5];
tm_ot=[210.04 216.5 ];

[n m]=size(jdy);
zsnk=ones(n,1)*.05;
ztsg=ones(n,1)*5;

%%%%%******  play with velocity structure functions

   cs=(cu+.75*cw)/2;%aver cu and cw for smoother structure function
   ii=find(isnan(cu));
	cs(ii)=.75*cw(ii);
	kk=find(isnan(cw));
   cs(kk)=cu(kk);
   
   
%**********************   scalar heat fluxes *******************************
wt=wt_son+0.51*(ta+tdk).*ustr0.*qstr0;  %correct sonic wt for humidity contribution 
hsc=wt*cpa.*rhoair; 
hsc_son=wt_son*cpa.*rhoair; 
hlc_lic=wq_lic*Le.*rhoair/1000;  
   

%*********  sub it tsg when snake is out of water
tsnk_m=tsg;zts=ztsg;
[n,m]=size(tm_in);
jg=[];
for k=1:m
   jf=find(jdy'>tm_in(k) & jdy'<tm_ot(k));%find when snake is in water
   jg=[jg jf];
end;
tsnk_m(jg)=tsnk(jg);%water temperature if snake is in
zts(jg)=zsnk(jg);%depth of water temp measurement
Ts=tsnk_m;
ii=find(isnan(tsnk));Ts(ii)=tsg(ii);
qse=qsea(Ts);

%**********************  get map of cruise track
lat=Lat;
lon=Lon;
%draw map of cruise track
  figure;map;
%**********************  get clear sky radiative fluxes
rlcler=5.67e-8*(ta+273.16).^4.*(.52+.085*sqrt(qa));%model clear sky
rlclr=(.52+.13/60*abs(lat)+(.082-.03/60*abs(lat)).*sqrt(qa)).*(5.67e-8*(ta+tdk).^4); 
rlcler=rlclr;
stdt=192;endt=floor(jdy(length(jdy)))+1;
% js=190;je=192;
iv=3;
AMMA_Rsclear_5;%get clear sky solar flux data


%**************   set instrument heights
zu=zu_etl(1);
zt=zt_etl(1);
zq=zq_etl(1);

jcool=1;
jwarm=1;
icount=1;
%*********************  housekeep variables  ********	
qcol_ac=0;
qjoule=0;
tau_ac=0;
jtime=0;
jamset=0;
tau_old=.06;
hs_old=10;
hl_old=100;
dsea=0;
dt_wrm=0;
tk_pwp=0;
fxp=.7;
ts_depth=.05;%bulk water temperature sensor depth

%******************  setup read data loop  **********	

[nn,nx]=size(x');  %detects # of lines of data  
ijx=1;
d1=datenum('31-Dec-2006');
clear dt;


for ibg = 1:nx 			%major read loop
	%***********   set variables not in data base  ********
	P=1010;                     %air pressure
	us=0;                       %surface current
	zi=600;                    %inversion ht
			%*******   decode date  ************
         
         
   [iyr mon iday ihr imin isec]=datevec(jdy(ibg)+d1);      
   
	%mon=fix(stux(1)/1e8)-iyr*100;
	%iday=fix(stux(1)/1e6)-iyr*1e4-mon*100;
	%ihr=fix(stux(1)/1e4)-iyr*1e6-mon*1e4-iday*100;
	%imin=fix(stux(1)/100)-fix(stux(1)/1e4)*100;
	%isec=0;
	timex=[iyr mon iday ihr imin isec];
	jd=jdy(ibg);
		%********   decode bulk met data ****
   u=U(ibg);%   
   phirel=reldir(ibg);%relative wind direction
	tsea=Ts(ibg);%bulk sea surface temp
	t=ta(ibg);%air temp
	qs=qse(ibg);%bulk sea surface humidity
	q=qa(ibg);%air humidity
	Rs=rs(ibg);%downward solar flux
	Rl=rl(ibg);%doward IR flux
   rain=org(ibg);%rain rate
   if sig_h<5
      jt=0;%flux data quality indicator
   else
      jt=1;
   end;

   tilt=tiltx(ibg);%flow tilt, fixed frame
   jm=sig_u(ibg);%ship manauver index
   lonn=Lon(ibg);%longitude
   latx=Lat(ibg);
  %  zu=17.5;%sensor height
  %  zq=zu;
  %  zt=zu;
   ts_depth=zts(ibg);%.05;%bulk water temperature sensor depth
  
   %*****  variables for warm layer  ***
	ntime=1e6*mon+1e4*iday+100*ihr+imin;
	time=(ihr*3600+imin*60)/24/3600;
	intime=time;
	loc=(lonn+7.5)/15;
	Rnl=.97*(5.67e-8*(tsea+273.16)^4-Rl);
	Rns=.945*Rs;
	%*********   set condition dependt stuff ******
	Le=(2.501-.00237*tsea)*1e6;
	cpv=cpa*(1+0.84*q/1000);
	rhoa=P*100/(Rgas*(t+tdk)*(1+0.61*q/1000));
	visa=1.325e-5*(1+6.542e-3*t+8.301e-6*t*t-4.8e-9*t*t*t);
	Al=2.1e-5*(tsea+3.2)^0.79;

	%**************   apply warm layer  ***********;
	if jwarm==1;                                            %do warm layer
		chktime=loc+intime*24;
		newtime=(chktime-24*fix(chktime/24))*3600;
		if icount>1                                     %not first time thru
			if newtime>21600 & jump==1
				%goto 16
			else
				if newtime <jtime		%re-zero at midnight
					jamset=0;
					fxp=.5;
					tk_pwp=19;
					tau_ac=0;
					qcol_ac=0;
					dt_wrm=0;
					jump=0;                   %goto 16
				else
						%****   set warm layer constants  ***
					rich=.65;                		%crit rich	
					ctd1=sqrt(2*rich*cpw/(Al*grav*rhow));
					ctd2=sqrt(2*Al*grav/(rich*rhow))/(cpw^1.5);
						%*********************************
					dtime=newtime-jtime;			%delta time for integrals
					qr_out=Rnl+hs_old+hl_old+RF_old;	%total cooling at surface
					q_pwp=fxp*Rns-qr_out;			%tot heat abs in warm layer
					if q_pwp<50 & jamset==0			%check for threshold
						%goto 16		
					else
						jamset=1;			%indicates threshold crossed
						tau_ac=tau_ac+tau_old*dtime;	%momentum integral
						if qcol_ac+q_pwp*dtime>0	%check threshold for warm layer existence
							for i=1:5		%loop 5 times for fxp
								fxp=1-(0.28*0.014*(1-exp(-tk_pwp/0.014))+0.25*0.357*(1-exp(-tk_pwp/0.357))+0.47*12.82*(1-exp(-tk_pwp/12.82)))/tk_pwp;
                        %fg=fpaul(tk_pwp);fxp=fg(1);
                        qjoule=(fxp*Rns-qr_out)*dtime;
								if qcol_ac+qjoule>0
									tk_pwp=min(19,ctd1*tau_ac/sqrt(qcol_ac+qjoule));
								end;
							end;%  end i loop
						else				%warm layer wiped out
							fxp=0.75;
							tk_pwp=19;
						end;%   end sign check on qcol_ac
						qcol_ac=qcol_ac+qjoule;		%heat integral
							%*******  compute dt_warm  ******
						if qcol_ac>0
							dt_wrm=ctd2*(qcol_ac)^1.5/tau_ac;
						else 
							dt_wrm=0;
						end;
					end;%                    end threshold check
				end;%                            end midnight reset
				if tk_pwp<ts_depth
					dsea=dt_wrm;
				else
					dsea=dt_wrm*ts_depth/tk_pwp;
				end;
						end;%                                    end 6am start first time thru
		else
			jump=1;					%idicates data start period
		end;%                                            end first time thru check
		jtime=newtime;
	end;%                                                     end warm layer model appl check
	
	
   ts=tsea+dsea;
   qs=qsea(ts);
	xx=[u us ts t qs q Rs Rl rain zi  P zu zt zq latx jcool 0 0 0] ;		%set data for basic flux alogithm
		%********    call modified LKB routine *******
	y=cor30a(xx);%latest version coare model
		%************* output from routine  *****************************
        hsb=y(1);                   %sensible heat flux W/m/m
        hlb=y(2);                   %latent
        tau=y(3);                   %stress
        zo=y(4);                    %vel roughness
        zot=y(5);                   %temp "
        zoq=y(6);                   %hum  "
        L=y(7);                     %Ob Length
        usr=y(8);                   %ustar
        tsr=y(9);                   %tstar
        qsr=y(10);                  %qstar  [g/g]
        dter=y(11);                 %cool skin delta T
        dqer=y(12);                 %  "   "     "   q
        tkt=y(13);                  %thickness of cool skin
        RF=y(14);                   %rain heat flux
        Wbar_webb=y(15);            %webb mean w 
        Cd=y(16);                   %drag @ zu
        Ch=y(17);                   %
        Ce=y(18);                   %Dalton
        Cdn_10=y(19);               %neutral drag @ 10 [includes gustiness]
        Chn_10=y(20);               %
        Cen_10=y(21);               %
        zax(1)=jd;                  %julian day
        zax(2:10)=xx(1:9);           %
        zax(4)=ts;                %Tsea [no cool skin]
        zax(11:32)=y;               %
        zax(33:35)=[dt_wrm tk_pwp ts];  %warm layer deltaT, thickness, corrected Tsea
	%*******   previous values from cwf hp basic code *****
	hsx=Hsb(ibg);
	hlx=Hlb(ibg);
   taux=Taub(ibg);
   
	%********************  save various parts of data **********************************
	%zax(35:38)=[hsx hlx taux zu];
   dt(ibg,:)=zax;
   %[jd  u  us  Tsea[no cool skin]  t qs q Rs  Rl rain hsb hlb tau zo zot zoq L  usr tsr qsr dter dqer tkt RF wbar Cd  Ch  Ce Cdn_10 Chn_10 Cen_10 ug dt_wrm tk_pwp  ts]
   %  1  2   3          4           5 6  7  8  9  10   11  12  13  14 15  16  17 18  19  20  21   22   23  24 25   26  27  28   29     30     31   32   33      34   35
   
   %struc(ibg,:)=[ct cq cu cw];
     
%   eyes(ibg,:)=[jt oph tilt jm phirel];
	e=[num2str(mon) '  ' num2str(iday) '  ' num2str(ihr) '   ' num2str(dter) '  ' num2str(dt_wrm) '  ' num2str(tk_pwp) ];
	disp(e)
	hs_old=hsb;
	hl_old=hlb;
	RF_old=RF;
	tau_old=tau;
   icount=icount+1;
   if phirel<90 | phirel>270  %save subset with good wind direction
      dtx(ijx,:)=zax;
      %idfx(ijx,:)=[hsix hlix tauix];   
      ijx=ijx+1;
   end;
end; %  data line loop

pst=psitg(zq./dt(:,17));

[np,a]=size(dt);
dtw=dt(:,33);
zw=dt(:,34);

Hsbb=dt(:,11);
Hlbb=dt(:,12);
HRF=dt(:,24);
Taubb=dt(:,13);


hlc=(hlc_lic);%
hl_webb=dt(:,25).*rhoair*Le/1000.*qa;


%****************Webb correction ozone fluxes*************
O3w_oc=O3w_o-dt(:,25);              % w'O3' calculated with crosscovariance, 10 min lagtime, webb correction
O3w_hc=O3w_h-dt(:,25);              % w'O3' calculated with crosscovariance, flux at hourly mean lagtime, webb correction
O3w_dc=O3w_d-dt(:,25);              % w'O3' calculated with crosscovariance, flux at dayly mean lagtime, webb correction
O3wc=O3w-dt(:,25);  

yfl(:,13)=O3w_oc;              % w'O3' calculated with crosscovariance, 10 min lagtime, webb correction
yfl(:,14)=O3w_hc;              % w'O3' calculated with crosscovariance, flux at hourly mean lagtime, webb correction
yfl(:,15)=O3w_dc;              % w'O3' calculated with crosscovariance, flux at dayly mean lagtime, webb correction
yfl(:,16)=O3wc;  
%*********************************************************

jdollars='AMMA';
 hl_lic=hlc_lic;
%create best covariance latent heat
hlc_bl=hlc;
% if ~isempty(hl_lic)
% 	ii=find(isnan(hlc));
% 	hlc_bl(ii)=hl_lic(ii);
% 	ii=find(isfinite(hlc) & isfinite(hl_lic));
% 	if ~isempty(ii);hlc_bl(ii)=(hlc(ii)+hl_lic(ii))/2;end;
% end;

AMMA_da_id2;
hlib_lic(ii1)=NaN;hlib_lic(ii2)=NaN;
%********************   ocean warm layer profiles  *******
dtz=zeros(np,6);
zz=0:5;
for i=1:6
   ii=find(zw>zz(i));
   dtz(ii,i)=dtw(ii).*(1-zz(i)./zw(ii));  %warming from mixed layer to z
end;
dt5=dtz(:,6);
tmxd=tsg-dtz(:,6);%mixed layer temp:  tsg with with 5-m warm layer removed

rnl=0.97*(5.67e-8*(Ts+273.15).^4-rl);
hnet=0.95*rs-rnl-Hsbb-Hlbb-HRF;
% stdt=192;
% endt=216;


%nar_diurnal;%perform diurnal averages, do plots

% stdt=192;
% endt=216;
AMMA_daily;%perform daily averages, do plots
 %prt_AMMA_flux_10;


AMMA_hr_ave;%convert to one hr averages
hr_bulk_AMMA;%recompute bulk, ID, and covariance fluxes at 1-hr
med_hl_plt_AMMA;
taucx=ytauc;
taucy=ytaucy;
tauib=ytauib;
taub=ytaub;
reldir=mod(thet1rel+360,360);

med_tau_plt_AMMA;
prt_AMMA_flux_hr;

%ii=find(lat>40 & isfinite(hsbb) & isfinite(hlbb));mean([ hsbb(ii) hlbb(ii) rs(ii) Rscl(ii) rl(ii) rlclr(ii) rnl(ii) hnet(ii)])