disp('bulk_tn')
%toga coare bulk flux model version 2.5
%***************************************
%uses following matlab subroutines:
%	julday.m
%	cool_lkb.m
%	psiu.m
%	psit.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 **************
%	YYMMHHMMSS=		date in toga coare 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)

%********** output data  ***************
%	hsb=			sensible heat flux (w/m^2)
%	hlb=			latent heat flux (w/m^2)
%	RF=			rain heat flux(w/m^2)
%	hl_webb=		webb addition to latent flux (w/m^2)
%	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');
flist=fopen('coartn_3.txt','r');
stuff=fgetl(flist);
stuff=fgetl(flist);

zu=15;
zt=15;
zq=15;
ts_depth=.05;

jcool=1;
jwarm=1;
icount=1;
%*********************  housekeep variables  ********	
qcol_ac=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;
%*******************  set constants  ****************
	tdk=273.16;
	grav=9.72;
	Rgas=287.1;
	cpa=1004.67;

	be=0.026;
	cpw=4000;
	rhow=1022;
	visw=1e-6;
	tcw=0.6;

%******************  setup read data loop  **********
nx=1620;				%# of lines of data     
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
	%*********************************************
	stuff=fgetl(flist);	    %read the data 1 line at a time
	%*********************************************
		%*******   decode date  ************
	iyr=eval(stuff(1:2));
	mon=eval(stuff(3:4));
	iday=eval(stuff(5:6));
	ihr=eval(stuff(7:8));
	imin=eval(stuff(9:10));
	isec=eval(stuff(11:12));
	timex=[iyr mon iday ihr imin isec];
	jd=julday(timex);
		%********   decode bulk met data ****
	u=eval(stuff(16:21));
	tsea=eval(stuff(38:42));
	t=eval(stuff(43:48));
	qs=eval(stuff(49:54));
	q=eval(stuff(55:60));
	Rs=eval(stuff(112:117));
	Rl=eval(stuff(118:123));
	rain=eval(stuff(124:129));
	lon=eval(stuff(212:220));
		%*****  variables for warm layer  ***
	ntime=1e6*mon+1e4*iday+100*ihr+imin;
	time=(ihr*3600+imin*60)/24/3600;
	intime=time;
	loc=(lon+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.27*0.357*(1-exp(-tk_pwp/0.357))+0.45*12.82*(1-exp(-tk_pwp/12.82)))/tk_pwp;
								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;
	x=[u us ts t qs q Rs Rl rain zi  P zu zt zq] ;		%set data for basic flux alogithm
		%********    call modified LKB routine *******
	y=cool_lkb(x);
		%************* 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
        hl_webb=y(15);              %webb contrib to latent heat flux
        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)=x(1:9);           %
        zax(4)=tsea;                %Tsea [no cool skin]
        zax(11:31)=y;               %
        zax(32:34)=[dt_wrm tk_pwp ts];  %warm layer deltaT, thickness, corrected Tsea
	%*******   previous values from cwf hp basic code *****
	hsx=eval(stuff(71:75));
	hlx=eval(stuff(86:91));
	taux=eval(stuff(105:111));
	%******************************************************
	zax(35:37)=[hsx hlx taux];
	dt(ibg,:)=zax;
	e=[stuff(16:21) stuff(38:60) stuff(112:123) stuff(71:75) stuff(86:90) stuff(105:111) '  ' num2str(hsb) '  ' num2str(hlb) '  ' num2str(tau) '  ' num2str(dt_wrm) '  ' num2str(tk_pwp)];% '  ' num2str(Cen_10)];
	disp(e)
	hs_old=hsb;
	hl_old=hlb;
	RF_old=RF;
	tau_old=tau;
	icount=icount+1;
end; %  data line loop
	
	
	
	
	
	
	
	
