disp('jasm_26')
%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');
%
didread=0;
if didread~=1;
   x=load('d:\cwf\flux_sum.txt');%read file with half-hour average data; set your local directory 
   jdy=x(:,1);%julian day at beginning of half hour average
   U=x(:,2); %true wind speed, m/s; etl sonic anemometer
   dir=x(:,3);%true wind direction (from), deg (clockwise rel north)
   tsnk=x(:,4);%sea snake temperature, C (0.05 m depth)
   tsg=x(:,5);%tsg water temperature, C (5 m depth)
   sal=x(:,6);%tsg salinity, psu (5 m depth)
   ta=x(:,7);%air temperature, C (z=14.5 m)
   qse=x(:,8);%sea surface specific humidity, g/kg (from snake when in water, see field 17)
   qa=x(:,9);%air specific humidity, g/kg (z=14.5  m)
   rs=x(:,10);%downward solar flux, W/m^2 (ETL units)
   rl=x(:,11);%downward IR flux, W/m^2 (ETL units)
   org=x(:,12);%rainrate, mm/hr (ETL STI optical rain gauge, uncorrected)
   ushp=x(:,13);%ship speed, m/s (SCS ODEC)
   head=x(:,14);%ship heading, deg clockwise rel north (SCS laser ring gyro)
   urel=x(:,15);%relative wind speed, m/s (ETL sonic)
   reldir=x(:,16);%relative wind direction (from), deg clockwise rel ship's bow(ETL sonic)
   zts=x(:,17);%depth for flux Ts reference, =0.05 when snake is in the water, 5.0 otherwise
   sig_u=x(:,18);%std dev of ship speed, m/s (>.1 indicates manauver during half hour average)
   Taub=x(:,19);%bulk wind stress along mean wind, N/m^2
   Hsb=x(:,20);%bulk sensible heat flux, W/m^2
   Hlb=x(:,21);%bulk latent heat flux, W/m^2 (includes Webb et al. correction)
   Hrain=x(:,22);%rain heat flux, W/m^2 as per Gosnell et al
   Lat=x(:,23);%latitude, deg  (SCS pcode)
   Lon=x(:,24);%longitude, deg (SCS pcode)
   rnl=x(:,25);%net longwave flux, W/m^2 using emissivity=0.97
   hnet=x(:,26);%net heat flux into ocean, W/m^2
   Ts=x(:,27);%Ts used in bulk flux calculation, =tsnk if zt=0.05
end;


zu=17.5;
zt=15.5;
zq=15.5;
ts_depth=.05;%bulk water temperature sensor depth

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=8080;				%# of lines of data  

[nn,nx]=size(x');
ijx=1;
d1=datenum('31-Dec-1998');

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);%wind speed
   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
   jt=0;%flux data quality indicator
   oph=9;%ophir data quality (bad=>15)
   tilt=5;%flow tilt, fixed frame
   jm=sig_u(ibg);%ship manauver index
   lon=Lon(ibg);%longitude
 %  zu=17.5;%sensor height
  % zq=zu;
  %zt=zu;
	ts_depth=zts(ibg);%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=(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.47*0.357*(1-exp(-tk_pwp/0.357))+0.25*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;
	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=cor26(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=Hsb(ibg);
	hlx=Hlb(ibg);
   taux=Taub(ibg);
   
	%********************  save various parts of data **********************************
	%zax(35:38)=[hsx hlx taux zu];
   dt(ibg,:)=zax;
   %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
[np,a]=size(dt);
dtw=dt(:,32);
zw=dt(:,33);
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



%********************  diurnal warm layer  ***********************
jj=find((jdy-floor(jdy))==0);
jj0=floor(jdy(1));
jjm=floor(jdy(nx));
tzil=tsg(jj);	
for mm=1:jjm-jj0+1
   im=(mm-1)*48+1;jm=im+47;
   t0(im:jm)=tzil(mm);
end;
figure;plot(jdy,tmxd,jdy,tsg);xlabel('Julian Day (1999)');ylabel('Temperature (C)');
Title('blu=mixed layer T, grn= 5m T (tsg)');

isn=find(zts==0.05);
figure;plot(jdy,dtw,jdy(isn),tsnk(isn)-tmxd(isn));xlabel('Julian Day (1999)');ylabel('DT Warm Layer (C)');
Title('blu=Coare26 Model, grn= Tsnake-tsg+DT5');


%****************  diurnal averages  *************
hr=0:23;hr=hr+.5;
dj=1/24;
jdx=jdy-floor(jdy);
for j=0:23
   ck=j*dj;
   i=j+1;
   m=find(jdx>ck & jdx<ck+dj);
   rainm(i)=mean(org(m));
   clr(i)=max(rs(m));
   rsm(i)=mean(rs(m));
   rlm(i)=mean(rl(m));
   rnlm(i)=mean(rnl(m));
	hnm(i)=mean(hnet(m));
   um(i)=mean(U(m));
   hlm(i)=mean(Hlb(m));
   hsm(i)=mean(Hsb(m));
   tsm(i)=mean(Ts(m));
   tsgm(i)=mean(tsg(m));
   tmxdm(i)=mean(tmxd(m));
   t5m(i)=mean(dt5(m));
   dtm(i)=mean(dtw(m));
end;

figure;plot(hr,rainm);xlabel('Hour of the Day');ylabel('Rainrate (mm/hr)');
figure;plot(hr,clr,hr,rsm,hr,hnm);xlabel('Hour of the Day');ylabel('Solar and Net Heat Flux (W/m^2)');
figure;plot(hr(1:12),rsm(1:12)./clr(1:12));xlabel('Hour of the Day');ylabel('Solar Cloud Transmission Coeff.');
figure;plot(hr,um);xlabel('Hour of the Day');ylabel('Wind Speed (m/s)');
figure;plot(hr,hlm,hr,-rnlm);xlabel('Hour of the Day');ylabel('Latent and Net IR Heat Fluxes (W/m^2)');
figure;plot(hr,rlm);xlabel('Hour of the Day');ylabel('Downward IR Flux (W/m^2)');
figure;plot(hr,tsm,hr,tsgm,hr,tmxdm+dtm,'.',hr,tsgm(1)+t5m,'x');xlabel('Hour of the Day');ylabel('Ocean Near Surface Temperatures (C)');

	
%************************  daily averages  ***********	
jj0=jdy(1);
[nn,nx]=size(jdy);
mx=floor(nn/48);
for i=1:mx
   ii=find(floor(jdy)==jj0);
   raind(i)=mean(org(ii));
	rsd(i)=mean(rs(ii));
   rld(i)=mean(rl(ii));
   rnld(i)=mean(rnl(ii));
	hnd(i)=mean(hnet(ii));
   ud(i)=mean(U(ii));
   hsd(i)=mean(Hsb(ii));
   hld(i)=mean(Hlb(ii));
   tsd(i)=mean(Ts(ii));
   tsgd(i)=mean(tsg(ii));
   tmxdd(i)=mean(tmxd(ii));
   t5d(i)=mean(dt5(ii));
   dtd(i)=mean(dtw(ii));
   jj0=jj0+1;   
end;

   
	
