disp('read_epicmonitor_10');
%This program reads the flux, ceilometer,and microwave radiometer data in
%and defines the variables 
%It also computes the IR and solar clear sky fluxes
%The 10h implies version h of the 10-min PACS data base (fall 1999 to fall 2004)

   y=load('c:\data\cwf\data\epic_01\metz_epic_10h.txt');%read file with 10-min average data; set your local directory 
   %x=y(1:20771,:);
   %x=y(1:24659,:);
   %x=y(1:30338,:);
   x=y(1:length(y),:);
   yr=x(:,1);%year of exp   
   jdy=x(:,2);%julian day at beginning of time average
   
   ushp=x(:,3);%doppler log, SCS (m/s)
   U=x(:,4);%true wind,ETL sonic (m/s)
   dir=x(:,5);%true wind direction, ETL sonic (deg)
   urel=x(:,6);%relative wind speed, ETL (m/s)
   reldir=x(:,7);%relative wind direction (from),clockwise rel ship's bow, ETL sonic (deg)
	head=x(:,8);%ship heading, deg clockwise rel north, SCS laser ring gyro (deg)
   tsnk=x(:,9);%sea snake temperature,ETL, 0.05 m depth (C)
	ta=x(:,10);%air temperature, ETL (C)   
   qse=x(:,11);%sea surface specific humidity,from snake (g/kg)
   qa=x(:,12);%air specific humidity, ETL (g/kg)
   hsc=x(:,13);%sensible heat flux, covariance, ETL sonic anemometer(W/m^2)
   hsib=x(:,14);%sensible heat flux, ID, ETL sonic anemometer(W/m^2)
   hsb=x(:,15);%bulk sensible heat flux, (W/m^2)
   hlc=x(:,16);%latent heat flux, covariance, (W/m^2)
   hlib=x(:,17);%latent heat flux, ID, (W/m^2)
   hlb=x(:,18);%bulk latent heat flux, W/m^2 (includes Webb et al. correction)
   taucx=x(:,19);%covariance streamwise stress, ETL sonic anemometer (N/m^2)
   taucy=x(:,20);%covariance cross-stream stress, ETL sonic anemometer (N/m^2)
   tauib=x(:,21);%ID streamwise stress, ETL sonic anemometer (N/m^2)
   taub=x(:,22);%bulk wind stress along mean wind, (N/m^2)
   rs=x(:,23);%downward solar flux, ETL units (W/m^2)
   rl=x(:,24);%downward IR flux,  ETL units (W/m^2)
   org=x(:,25);%rainrate, ETL STI optical rain gauge, uncorrected (mm/hr)
   J=x(:,26);%ship plume contamination index
   sigoph=x(:,27);%standard deviation of ophir fast hygrometer clear channel
   tiltx=x(:,28);%flow tilt at ETL sonic anemometer, earth frame
   Jm=x(:,29);%ship maneuver index
   ct=x(:,30);%ct^2 (K^2/m^.667)
   cq=x(:,31);%cq^2 ((g/kg)^2/m^.667) 
	cu=x(:,32);%cu^2 ((m/s)^2/m^.667)
   cw=x(:,33);%cw^2 ((m/s)^2/m^.667)
   hrain=x(:,34);%rain heat flux, as per Gosnell et al 1995, JGR, 18437-18442, (W/m^2)
   hlwebb=x(:,35);%correction to measured latent heat flux, Webb et al. 1980,QJRMS, 85-100
   lat=x(:,36);%latitude, deg  (SCS pcode)
   lon=x(:,37);%longitude, deg (SCS pcode)
	zu_etl=x(:,38);%height of mean wind sensor, 17.7 m
   zt_etl=x(:,39);%height of mean air temperature sensor, 15.5 m
   zq_etl=x(:,40);%height of mean air humidity sensor, 15.5 m
   %*****   ships imet and scs data
   sog=x(:,41);%speed over ground, SCS gps, (m/s)
   U_scs=x(:,42); %true wind speed, imet propvane anemometer (m/s)
   dir_scs=x(:,43);%true wind direction (from),clockwise rel north, imet,(deg)
   cog=x(:,44);%%course over ground, SCS gps, (m/s)
   tsg=x(:,45);%tsg water temperature, 5 m depth, (C)
   ta_im=x(:,46);%imet air temperature (C) 
   qs_tsg=x(:,47);%imet bulk water specific humidity (g/kg) 
   qa_im=x(:,48);%imet air specific humidity, (g/kg)
   jdq=jdy+datenum(yr+1900,1,1)-datenum(1999,1,1);%(yr-99)*366+jdy;%set single JD base 1999

   %*************   correct some data problems
ii=1892;jj=2655;
qse(ii:jj)=qsea(tsnk(ii:jj));   
   
ib=find(tsnk<10);tsnk(ib)=tsg(ib); qse(ib)=qsea(tsnk(ib));
ib=find(tsnk>35);tsnk(ib)=tsg(ib); qse(ib)=qsea(tsnk(ib));
jj=find(rs<-100);rs(jj)=NaN;

%fix bad true wind dir in spring 02 data
rdcon=pi/180;
ode=ushp;
lrgcm=cos(head*pi/180);
lrgsm=sin(head*pi/180);
uu=-urel.*cos(reldir*pi/180)/.95;
vv=urel.*sin(reldir*pi/180)/1.15;
sodx=ode.*lrgcm;
sody=ode.*lrgsm;
un=uu.*lrgcm+vv.*lrgsm;
ue=uu.*lrgsm-vv.*lrgcm;
un1=un+sodx;
ue1=ue+sody;
dirx=mod(atan2(-(ue+sody),-(un+sodx))/rdcon+360,360);
ll=find(jdq>1100);
dir(ll)=dirx(ll);

%**************************

%**************   read ceilo data
clear xh yh yrc jdyc allda nocld onecld mltcld obscure pobs clrf ceilof ceilofobs zb zb1 zb2
yh=load('c:\data\cwf\data\epic_01\metz_epic_ceilo_10h.txt');%read file with 10-min average data; set your local directory 
   yrc=yh(:,1);%year of exp   
   jdyc=yh(:,2);%julian day at beginning of time average
   jdqc=datenum(yrc,1,1)-datenum(1999,1,1)+jdyc;%(yrc-1999)*366+jdyc;
   xh=interp1(jdqc,yh,jdq);
   allda=xh(:,3);%number of data points
   nocld=xh(:,4);%numbewr of clear
   onecld=xh(:,5);%number with one cloud
   mltcld=xh(:,6);%number with multiple clouds
   obscure=xh(:,7);%number obscured cases
	pobs=xh(:,8);%number partially obscured
   clrf=xh(:,9);%clear fraction
	ceilof=xh(:,10);%cloud fraction
   ceilofobs=xh(:,11);%cloud fraction, including partially obscure
   zb=xh(:,12);%median cloud base height (m)
   zb1=xh(:,13);%15% cloud base height (m)
   zb2=xh(:,14);%85% cloud base height (m)
   

%**************   read microwave data
clear xi yi yrm jdym vap liq xnm sigvap sigliq
yi=load('c:\data\cwf\data\epic_01\metz_epic_micro_10h.txt');%read file with 10-min average data; set your local directory 
   
yrm=yi(:,1);%year of exp   
   jdym=yi(:,2);%julian day at beginning of time average
   jdqm=datenum(yrm,1,1)-datenum(1999,1,1)+jdym;%(yrm-1999)*366+jdym;
   ibx=find(yi(:,4)>.6 & jdqm>920 & jdqm<1030);yi(ibx,4)=NaN;
   ibx=find(yi(:,7)>.6 & jdqm>920 & jdqm<1030);yi(ibx,7)=NaN;
   imx=find(~isnan(jdqm) & yi(:,6)<.4);
   xi=interp1q(jdqm(imx),yi(imx,:),jdq);
   vap=max(1,xi(:,3));%precipitable water (cm)
   liq=xi(:,4);%column liquid water (cm)
   xnm=xi(:,5);%
   sigvap=xi(:,6);%stand deviation, vap (cm)
   sigliq=xi(:,7);%stand deviation, liq (cm)
	
   %*****************
     
%*******************  set constants  ****************
tdk=273.16;
grav=grv(lat);%9.82;
Rgas=287.1;
cpa=1004.67;
rhoair=1010*100./(Rgas*(ta+tdk).*(1+0.61*qa/1000));%density of air, kg/m^3

ii=find(isfinite(tsnk));Le=(2.501-.00237*median(tsnk(ii)))*1e6;
rlclr=(.5+.13/60*abs(lat)+(.091-.03/60*abs(lat)).*sqrt(qa)).*(5.67e-8*(ta+tdk).^4); %2-param clear model
rlclr3=(.5+.13/60*abs(lat)+(.091-.03/60*abs(lat)).*sqrt(qa) -.0188+.0063*vap).*(5.67e-8*(ta+tdk).^4); %3-param clear model

js=320;%pick one day for solar plot

%fall solar parameters
p=1010;%pressure
k1=.001;%aerosol optical depth, band 1
k2=.001;%aerosol optical depth, band 2
oz=0.2;
ivz=qa/3.8;%fall b value
epic_Rsclear_10e;%get clear sky solar flux data

%identify cruise sections
ii=1:3817;%fall 99
jj=3818:7573;%spring 00
kk=7574:11288;%fall 00
ll=11299:15875;%spring 01
mm=15876:26450;%fall 01
nn=26451:30338;%spr 02
pp=30339:34082;%fall 02
pq=34083:36926;%fall 03
pr=36927:41210;%fall 04

rsclr=Rscl;
rsclr0=Rscl;%no aerosol

%****************   Adjust aerosol depth for each cruise

k1=.07;%aerosol optical depth, band 1
k2=.07;%aerosol optical depth, band 2
epic_Rsclear_10e
rsclr(pp)=Rscl(pp);
rsclr(pq)=Rscl(pq);

k1=.04;%aerosol optical depth, band 1
k2=.04;%aerosol optical depth, band 2
epic_Rsclear_10e
rsclr(pr)=Rscl(pr);


%spring solar parameters
p=1010;%pressure
k1=.04;%aerosol optical depth, band 1
k2=.04;%aerosol optical depth, band 2
oz=0.2;
latz=min(lat,12);latz=max(latz,-10);
ivz=3.5+2*(10-latz)/20;%spring b value
epic_Rsclear_10e
%insert spring values
rsclr(jj)=Rscl(jj);

k1=.09;%aerosol optical depth, band 1
k2=.09;%aerosol optical depth, band 2
epic_Rsclear_10e
rsclr(ll)=Rscl(ll);

k1=.07;%aerosol optical depth, band 1
k2=.07;%aerosol optical depth, band 2
epic_Rsclear_10e
rsclr(nn)=Rscl(nn);
%****************************************


rlnet=0.97*(rl-5.67e-8*(tsnk+tdk).^4);%net IR flux
rlnetc=0.97*(rlclr-5.67e-8*(tsnk+tdk).^4);%net IR clear sky flux
hnet=0.95*rs-rlnet-hsb-hlb-hrain;% net heat into the ocean

%find clear periods with good data, in specified latitude band, skipping
%days 661-664
ii=find(~isnan(rl(1:length(ceilof)-2)) & ~isnan(vap(1:length(ceilof)-2)) & (jdqc(1:length(ceilof)-2)<661 | jdqc(1:length(ceilof)-2)>664) & ceilof(1:length(ceilof)-2)==0 & abs(lat(1:length(ceilof)-2))<15);
jj=find(ceilof(ii-1)==0 & ceilof(ii+1)==0 & ceilof(ii-2)==0 & ceilof(ii+2)==0);length(jj)
%plot clear sky model vs observed flux in 'clear' periods
figure;h=plot(rlclr3(ii(jj)),rl(ii(jj)),'.',[320 440],[320 440]);set(gca,'fontsize',14);set(h,'LineWidth',2);axis([320 440 320 440]);xlabel('Clear Model IR Flux (Wm^{-2})');ylabel('Clear Measured IR Flux (Wm^{-2})');
figure;h=plot(rsclr(ii(jj)),rs(ii(jj)),'.',[0 1200],[0 1200]);set(gca,'fontsize',14);set(h,'LineWidth',2);axis([0 1200 0 1200]);xlabel('Clear Model Solar Flux (Wm^{-2})');ylabel('Clear Measured Solar Flux (Wm^{-2})');

figure;h=plot(U,taub,'.');set(gca,'fontsize',14);set(h,'LineWidth',2);axis([0 15 0 0.5]);xlabel('Wind Speed (m/s)');ylabel('Bulk Turbulent Stress (Nm^{-2})');
figure;h=plot(U,hlb,'.');set(gca,'fontsize',14);set(h,'LineWidth',2);axis([0 15 0 250]);xlabel('Wind Speed (m/s)');ylabel('Bulk Turbulent Latent Heat Flux (Wm^{-2})');


cloudliq_biascorrect_10h;
liq_nb=liq-liq0;%remove bias from cloud lwc
