% qflux.m - routine to generate heat flux for model runs

% History:
%  
%  first written & used at NCSC OSSE camp on EN302 cruise
%  second used by Jmanning on OC303 cruise

% WARNING:
%  Each cruise may have different input data so they
%  are taken individually. The cruises are dealt with
%  separately and are listed cronologically below
 

wc=menu('Which Cruise?','June 1996',...
                        'May 1997 (OC303)',...
                        'June 1997 (EN302)',...
                        'Feb 1997 (OC298)',...
                        'Feb 1999 (OC336)');

if wc==1;
  % June 1996 (using beardsley calculations from mooring)
  load sf4met1h.mat
  figure(1)
  plot(YD-1,Qnet)
  axis([160 172 -200 800])
  hold on;
  % filter radiation flux
  fqswlw=plfilt(Qsw+Qlw);
  plot(YD-1,fqswlw,'m','linewidth',3) 
  plot(YD-1,Qnlp,'g','linewidth',3)
  ylabel('W/m2 Net Heat Flux')
  plot(YD-1,Qlat+Qsen,'c','linewidth',3)
  legend('Qnet','Qsw+Qlw','Qnlp','Qlat+Qsen')
  ydcalab(96);


elseif wc==2;
  %  May 1997 (oceanus ship data and daily paine solor flux)
  load /home3/ocn/jmanning/gbs97/hcast/fcast/oc_shipdata
  % subsample, if necessary, and calculate wind_spd
  ssr=6;
  tl=length(ewind);
  k=0;
  disp('Caculating wind spd...');
  for i=1:ssr:tl;
    k=k+1;
    [spd(k),dir(k)]=uv2sd(ewind(i),nwind(i));
  end;

  disp('Warning: madeup constant humidity = 90')
  humidity=ones(length(spd),1)*90;

  % calculate Qnet=Qlat+Qsen+Qlw+Qsw using Bob B's "Air_Sea"
  % HFBULKTC,SWHF,LWHF
  path(path,'/home3/ocn/jmanning/gbs99/rtda/qflux')
  disp('Calculating Latent & Sensible...');
  A=hfbulktc(spd,10,temp_air(1:ssr:tl),10,humidity,10,...
     press_bar(1:ssr:tl),temp_ss5(1:ssr:tl));
   save qlatsen97 yrday0_gmt A;
  % load qlatsen97
   plot(yrday0_gmt(1:ssr:tl),A(:,1)+A(:,2),'m');
   hold on;

  % add daily solar flux of Payne 
  load /data5/jmanning/heatflux/payne567_v4
  ydd=dayno(732:length(dayno))-1.0;
  qsol=soradm(732:length(dayno))/24;
  plot(ydd,qsol,'r--','linewidth',2)

  % add hourly heat flux
  load /data5/jmanning/heatflux/payne7hr.dat
  phr=payne7hr;
  % note I added 4/24 to get gmt and -1 to get yrday0_gmt
  yd=(phr(:,1)+phr(:,2)/2400)-1+4/24;
  plot(yd,phr(:,3),'g--');

  % calculate hourly average from OCEANUS and add to Payne phr(:,3)
  disp('Calculating hourly average OCEANUS data ...');
  [yf,ff]=runave(yrday0_gmt(2:ssr:tl),A(:,1)+A(:,2),60,0,0,1);
  save yfff yf ff;
  %load yfff
  % interpolate this on to the hourly figure
  idyd=find(yd>126.2&yd<142.5);
  syf=interp1(ff,yf,yd(idyd));
  plot(yd(idyd),syf+phr(idyd,3),'b','linewidth',3);
  axis([126 132 -400 1200])
  ylabel('Heat Flux (W/m^2)')
  yd2calab(97)
  title('OC303 SHIP DATA & Woods Hole Insolation','fontsize',18);
  legend('Latent & Sensible',...
      'Daily Averaged Solar Flux',...
      'Hourly Averaged Solar Flux',...
      'hourly net')

  % save output
  id=find(yd(idyd)>=126&yd(idyd)<132);
  yy=[yd(idyd(id)) syf(id)+phr(idyd(id),3)];
  save /data5/jmanning/heatflux/oc303.dat yy -ascii

elseif wc==3;
  %  June 1997 (endeavor ship & payne daily flux)
  load /data5/jmanning/ship/en302/emet_W1_1997
  figure(2);
  % calculate Qnet=Qlat+Qsen+Qlw+Qsw using Bob B's "Air_Sea"
  % HFBULKTC,SWHF,LWHF
  %path(path,'bobb')
  %A=hfbulktc(wind_speed_c,10,temp_air,10,humidity,10,press_bar,temp_ss1)
  %save qlatsen97 yrday_utc A;
  load qlatsen97
  plot(yrday_utc-1,A(:,1)+A(:,2),'m');
  hold on;

  % add daily solar flux of Payne for plotting only
  load payne567
  ydd=dayno(732:length(dayno))-1.0;
  qsol=soradm(732:length(dayno))/24;
  plot(ydd,qsol,'r--','linewidth',2)

  % add hourly heat flux
  load payne7hr.dat
  phr=payne7hr;
  % note I added 4/24 to get gmt and -1 to get yrday0_gmt
  yd=(phr(:,1)+phr(:,2)/2400)-1+4/24;
  plot(yd,phr(:,3),'g--');

  % calculate hourly average from endeavor and add to Payne phr(:,3)
  disp('Calculating hourly average endeavor data ...');
  %[yf,ff]=runave(yrday_utc-1,A(:,1)+A(:,2),60,0,0,1);
  %save yfff yf ff;
  load yfff
  % interpolate this on to the hourly figure
  syf=interp1(ff,yf,yd);
  %plot(yd,syf,'m','linewidth',3);

  plot(yd,syf+ phr(:,3),'b','linewidth',3);
  axis([160 172 -100 1200])
  ylabel('Heat Flux (W/m^2)')
  yd2calab(97)
  title('EN302 Alongtrack Data & Woods Hole Insolation','fontsize',18);
  legend('Latent & Sensible',...
      'Daily Averaged Solar Flux',...
      'Hourly Averaged Solar Flux',...
      'hourly net')

  % save output
  id=find(yd>=160&yd<173);
  yy=[yd(id) syf(id)+phr(id,3)];
  %save /globec/manning/met/heatflux.dat yy -ascii

elseif wc==4;
  %  OC298 case (no ship data available) Feb 1997
  %  Load buoy 11 data
  load /data5/jmanning/wind/mtl/b1197.mtl
  d=b1197;

  % filter out bad data now
  id=find(d(:,6)~=99.0&d(:,13)~=999.0&d(:,14)~=999.0&d(:,12)~=9999.0);
  d=d(id,:);

  % calculate yearday
  for i=1:length(d);
    yrday0_gmt(i)=julian(d(i,2),d(i,3),d(i,1),d(i,4)*100);
  end;

  spd=d(:,6);
  temp_air=d(:,13);
  temp_ss5=d(:,14);
  press_bar=d(:,12);
  
  %subsample rate
  ssr=6;
  tl=length(d);

  disp('Warning: madeup constant humidity = 90')
  humidity=ones(length(spd),1)*90;

  % calculate Qnet=Qlat+Qsen+Qlw+Qsw using Bob B's "Air_Sea"
  % HFBULKTC,SWHF,LWHF
  path(path,'/home3/ocn/jmanning/atsee/qflux')
  %plot input variables 
  %  pltinpvar;
  disp('Calculating Latent & Sensible...');
  A=hfbulktc(spd(1:ssr:tl),10,temp_air(1:ssr:tl),10,...
     humidity(1:ssr:tl),10,...
     press_bar(1:ssr:tl),temp_ss5(1:ssr:tl));
   save qlatsenF97 yrday0_gmt A;
  % load qlatsen97
   plot(yrday0_gmt(1:ssr:tl),A(:,1)+A(:,2),'m');
   hold on;

  % add daily solar flux of Payne 
  load /data5/jmanning/heatflux/payne567_v4
  ydd=dayno(732:length(dayno))-1.0;
  qsol=soradm(732:length(dayno))/24;
  plot(ydd,qsol,'r--','linewidth',2)

  % add hourly heat flux (NONE AVAILABLE)
  %load /data5/jmanning/heatflux/payne7hr.dat
  %phr=payne7hr;
  % note I added 4/24 to get gmt and -1 to get yrday0_gmt
  %yd=(phr(:,1)+phr(:,2)/2400)-1+4/24;
  %plot(yd,phr(:,3),'g--');
  yd=31:1/24:58;

  % calculate hourly average from OCEANUS and add to Payne phr(:,3)
  disp('Calculating hourly average OCEANUS data ...');
  [yf,ff]=runave(yrday0_gmt(2:ssr:tl),A(:,1)+A(:,2),60,0,0,1);
  save yfff yf ff;
  %load yfff
  % interpolate this on to the hourly figure
  idyd=find(yd>31&yd<58);
  syf=interp1(ff,yf,yd(idyd));
%  plot(yd(idyd),syf+phr(idyd,3),'b','linewidth',3);
  plot(yd(idyd),syf,'b','linewidth',3);
  axis([31 58 -200 200])
  ylabel('Heat Flux (W/m^2)')
  yd2calab(97)
  title('OC298 "ESTIMATE"','fontsize',18);
  text(min(yd),150,'using NOAA BUOY 44011 & Woods Hole Insolation','fontsize',14);
  legend('Latent & Sensible',...
      'Daily Averaged Solar Flux',...
      'net')
%      'Hourly Averaged Solar Flux',...

  % save output
  id=find(yd(idyd)>=31&yd(idyd)<58);
  yy=[yd(idyd(id))' syf(id)];
  save /data5/jmanning/heatflux/oc298.dat yy -ascii

elseif wc==5;
  %  OC336 case ( ship data available) Feb 1999
  %  Load  data
  load /data5/jmanning/ship/oc336/oc336.mtl
  d=oc336;

  % filter out bad data now
  %id=find(d(:,6)~=99.0&d(:,13)~=999.0&d(:,14)~=999.0&d(:,12)~=9999.0);
  %d=d(id,:);

  % calculate yearday
  %for i=1:length(d);
  %  yrday0_gmt(i)=julian(d(i,2),d(i,3),d(i,1),d(i,4)*100);
  %end;
  yrday0_gmt=d(:,1)-1;

  spd=d(:,13);
  temp_air=d(:,9);
  temp_ss5=d(:,5);
  press_bar=d(:,7);
  
  %subsample rate
  ssr=6;
  tl=length(d);

  %disp('Warning: madeup constant humidity = 90')
  %humidity=ones(length(spd),1)*90;
  humidity=d(:,10);

  % calculate Qnet=Qlat+Qsen+Qlw+Qsw using Bob B's "Air_Sea"
  % HFBULKTC,SWHF,LWHF
  path(path,'/home3/ocn/jmanning/atsee/qflux')
  %plot input variables 
  %  pltinpvar;
  disp('Calculating Latent & Sensible...');
  A=hfbulktc(spd(1:ssr:tl),10,temp_air(1:ssr:tl),10,...
     humidity(1:ssr:tl),10,...
     press_bar(1:ssr:tl),temp_ss5(1:ssr:tl));
   save qlatsenF99 yrday0_gmt A;
  % load qlatsen99
   plot(yrday0_gmt(1:ssr:tl),A(:,1)+A(:,2),'m');
   hold on;

  % add daily solar flux of Payne 
  %load /data5/jmanning/heatflux/payne567_v4
  %ydd=dayno(732:length(dayno))-1.0;
  %qsol=soradm(732:length(dayno))/24;
  qsol=d(:,8);
  %  plot(ydd,qsol,'r--','linewidth',2)
    plot(yrday0_gmt,qsol,'r--','linewidth',2)

  % add hourly heat flux (NONE AVAILABLE)
  %load /data5/jmanning/heatflux/payne7hr.dat
  %phr=payne7hr;
  % note I added 4/24 to get gmt and -1 to get yrday0_gmt
  %yd=(phr(:,1)+phr(:,2)/2400)-1+4/24;
  %plot(yd,phr(:,3),'g--');
  %yd=40:1/24:53;
  % calculate hourly average from OCEANUS and add to Payne phr(:,3)
  %disp('Calculating hourly average OCEANUS data ...');
  %[yf,ff]=runave(yrday0_gmt(2:ssr:tl),A(:,1)+A(:,2),60,0,0,1);
  %save yfff yf ff;
  %load yfff
  % interpolate this on to the hourly figure
  %idyd=find(yd>40&yd<53);
  %syf=interp1(ff,yf,yd(idyd));
  %  plot(yd(idyd),syf+phr(idyd,3),'b','linewidth',3);
  %plot(yd(idyd),syf,'b','linewidth',3);

  % add to get net
  net=A(:,1)+A(:,2)+qsol(1:ssr:length(qsol));
  plot(yrday0_gmt(1:ssr:length(d)),net,'b','linewidth',3);
  axis([40 53 -400 600])

  ylabel('Heat Flux (W/m^2)')
  yd2calab(99)
  title('OC336 "ESTIMATE"','fontsize',18);
  text(min(yd),150,'using OCEANUS alongtrack data','fontsize',14);
  legend('Latent & Sensible',...
      'Short wave Radiation',...
      'net')
%      'Hourly Averaged Solar Flux',...

  % save output
  %id=find(yd(idyd)>=40&yd(idyd)<43);
  yy=[yrday0_gmt(1:ssr:length(d)) net];
  save /data5/jmanning/heatflux/oc336.dat yy -ascii

end;