% Program written to incorporate CF warm layer code in coare35vn bulk
% algorithm
% Ludovic Bariteau, Jim Edson & Chris Fairall
% 1st version: March 2013

clear all;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grab the 10 minute file from DYNAMO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load Revelle10minutesLeg3_r3; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs Required
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Yday           %decimal yearday
%Ur             %relative wind speed, m/s
%Tsea           %Sea temperature, C at ts_depth 
%Tair           %air temperature, C
%Pair           %air temperature, C
%RH             %relative humidity (%)
%Lat            %latitude, deg
%Lon;           %longitude, deg
%Solar          %downwelling solar flux, W/m^2 (negative down)
%IR             %downwelling IR flux, W/m^2    (negative down)
%Rainrate       %rainrate, mm/hr

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assign variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zu=10;              %anemometer height
zt=10;              %air T height
zq=10;              %humidity height

Yday=yday;          %decimal yearday Jan 1 at noon = 1.5
Ur=Ur10;            %relative wind speed
Pair=Pair10;        %pressure, mb
Tair=T10;           %air temperature, C
RH=RH10;            %relative humidity, %
Tsea=Tsea-dsea;     %bulk sea temperature, C - The saved value has already
                    %been adjusted for the warm layer.  Removal of dsea
                    %provides the measured value, which should be passed 
                    %to the warm layer code.
dseachk=dsea;
ts_depth=0.05;      %bulk sea temperature sensor depth
Solar=-Solardn;     %downward solar flux, W/m^2
IR=-IRdn;           %downward IR flux, W/m^2
Rainrate=P;         %rainrate, mm/hr
Lat=Lat;            %latitude, deg
Lon=Lon;            %longitude, deg
zi=600;             %inversion ht

%************************************************************
%  This test program uses the sea snake temperature at
%  approximately 5 cm to compute the depth and amplitude of 
%  the warm layer.  The heating at this depth is added to 
%  sea snake and the cool-skin affect is removed in the 
%  coare35vn algorithm.  Therefore, fluxes in the returned 
%  B array have been computed based on a Tsea measurements 
%  that have been corrected for the warm layer and cool skin.
%************************************************************

B=coare35vnWarm(yday,Ur,zu,Tair,zt,RH,zq,Pair,Tsea,Solar,IR,Lat,Lon,zi,Rainrate,ts_depth);

%Outputs
%B(1)-B(39):     %See COARE 3.5 output list
dt_wrm=B(:,40)'; %warming (deg C) across entire warm layer
dz_wrm=B(:,41)'; %warm layer thickness (m)
dsea=B(:,42)';   %heating (deg C) at selected depth (i.e., ts_depth)

%************************************************************
% Make sure values are consistent
%************************************************************

figure(1);clf
plot(dsea,dseachk,'.',[0 0.15],[0 0.15],'k')
xlabel('\Delta T new (^oC)','fontsize',18)
ylabel('\Delta T old (^oC)','fontsize',18)

%************************************************************
%  The following shows how to use the warm layer variables to
%  compare the temperature measured at another depth, i.e.,
%  the ThermoSaliniGraph (TSG) on the Revellle. 
%  If the sensor is beneath the warm layer depth, the total 
%  temperature change (the warm layer amplitude) is added to 
%  the sensor. 
%  Note that these are the values just below the surface 
%  without the cool skin correction.
%************************************************************

tsg_depth=4;                      %Approximate depth of TSG
dseag=dt_wrm./dz_wrm.*tsg_depth;  %Linear heating profile is assumed
                                  % dT(z) = dT/dz*z
i=find(dz_wrm<tsg_depth);         %Use the total warming if the sensor
dseag(i)=dt_wrm(i);               %is below the warm layer.

tsg=TseaTSG+dseag;                %Add warm layer correction to adjust
ts=Tsea+dsea;                     %temps to value just below surface.

%************************************************************
% The expectation is that the TSG values will give
% better agreement with the Sea Snake after the warm layer
% correction has been applied.
%************************************************************

figure(2);clf
plot(yday,Tsea,'m',yday,ts,'r','linewidth',2)
xlabel('Yearday 2011','fontsize',18)
ylabel('Temperatures (^oC)','fontsize',18)
hold
plot(yday,TseaTSG,'b',yday,tsg,'k','linewidth',2)
axis([312 342 28.5 33])
legend('Sea Snake at 0.05 m','Sea Snake Corrected to surface','TSG at ~4.0 m','TSG Corrected to surface',1)
hold

