%This script computes both sides of the thermal wind relation using NCEP
%Reanalysis climatological data.
%------------------------------------------------------
%USER DEFINED INPUT VARIABLES

%INPUT DATA FILE NAMES
fname1='t_500.nc';
fname2='u_400.nc'; %For the wind, put the lower pressure level first
fname3='u_600.nc';

%INPUT PRESSURE LEVEL (mb) FOR TEMPERATURE CROSS-SECTION

p=500;

%INPUT PRESSURE DIFFERENCE (mb) BETWEEN TWO WIND LEVELS

dp=-200;

%INPUT SOUTHERN LATITUDE FOR CROSS-SECTION
%Note: This script is for degrees North latitude in 2.5 degree increments from 0 to 90N.
%Note: Cross-section by default will cut through 25 degrees of latitude, but you can modify the
%script to increase or decrease the cross-section length.
%Note: Southern Hemisphere may need latitude adjustment.  Plot the data to
%see.

lat=25;

%INPUT LONGITUDE FOR CROSS-SECTION
%Note: Please write a number divisible by 2.5, with 0 the Prime Meridian
%and 357.5 the equivalent of 2.5 W.

lon=150;

%------------------------------------------------------

%Read NetCDF files

nc=netcdf(fname1,'nowrite');
%ncdump(nc);

nc2=netcdf(fname2,'nowrite');
%ncdump(nc2);

nc3=netcdf(fname3,'nowrite');
%ncdump(nc3);

%------------------------
%Now, we need to specify the longitude, and latitude range.  The below
%commands are in (time,level,lat,lon) format, with a maximum value of (1 1
%37 144).  That is, 1 time step (your chosen month), 1 level (your pressure
%level), 37 possible latitudes at 2.5 degree spacing (1=90N, 37=0N), and
%144 possible longitudes also at 2.5 degree spacing (1=0E, 144=2.5W).

%The below example takes a cross-section at 150E (lon=61) from 25-50N (lat=27 to lat=17).
x=((90-lat)/2.5)+1;
y=(lon/2.5)+1;
R=287.06; %Units are J/kg/K
f=2*0.0000729*sin(lat*pi/180); %Units are in 1/s

temp0=nc{'air'}(1,1,x+1,y);
temp=nc{'air'}(1,1,x,y);
temp2=nc{'air'}(1,1,x-1,y);
temp3=nc{'air'}(1,1,x-2,y);
temp4=nc{'air'}(1,1,x-3,y);
temp5=nc{'air'}(1,1,x-4,y);
temp6=nc{'air'}(1,1,x-5,y);
temp7=nc{'air'}(1,1,x-6,y);
temp8=nc{'air'}(1,1,x-7,y);
temp9=nc{'air'}(1,1,x-8,y);
temp10=nc{'air'}(1,1,x-9,y);
temp11=nc{'air'}(1,1,x-10,y);
temp12=nc{'air'}(1,1,x-11,y);

xrange=[lat lat+2.5 lat+5 lat+7.5 lat+10 lat+12.5 lat+15 lat+17.5 lat+20 lat+22.5 lat+25];

temp_all=[temp temp2 temp3 temp4 temp5 temp6 temp7 temp8 temp9 temp10 temp11];
temp_grad=(1/555000)*[temp2-temp0 temp3-temp temp4-temp2 temp5-temp3 temp6-temp4 temp7-temp5 temp8-temp6 temp9-temp7 temp10-temp8 temp11-temp9 temp12-temp10];
tw_rhs=R*temp_grad/f/p;


uwind=nc2{'uwnd'}(1,1,x,y);
uwind2=nc2{'uwnd'}(1,1,x-1,y);
uwind3=nc2{'uwnd'}(1,1,x-2,y);
uwind4=nc2{'uwnd'}(1,1,x-3,y);
uwind5=nc2{'uwnd'}(1,1,x-4,y);
uwind6=nc2{'uwnd'}(1,1,x-5,y);
uwind7=nc2{'uwnd'}(1,1,x-6,y);
uwind8=nc2{'uwnd'}(1,1,x-7,y);
uwind9=nc2{'uwnd'}(1,1,x-8,y);
uwind10=nc2{'uwnd'}(1,1,x-9,y);
uwind11=nc2{'uwnd'}(1,1,x-10,y);


uwind_2=nc3{'uwnd'}(1,1,x,y);
uwind2_2=nc3{'uwnd'}(1,1,x-1,y);
uwind3_2=nc3{'uwnd'}(1,1,x-2,y);
uwind4_2=nc3{'uwnd'}(1,1,x-3,y);
uwind5_2=nc3{'uwnd'}(1,1,x-4,y);
uwind6_2=nc3{'uwnd'}(1,1,x-5,y);
uwind7_2=nc3{'uwnd'}(1,1,x-6,y);
uwind8_2=nc3{'uwnd'}(1,1,x-7,y);
uwind9_2=nc3{'uwnd'}(1,1,x-8,y);
uwind10_2=nc3{'uwnd'}(1,1,x-9,y);
uwind11_2=nc3{'uwnd'}(1,1,x-10,y);

tw_lhs=(1/dp)*[uwind-uwind_2 uwind2-uwind2_2 uwind3-uwind3_2 uwind4-uwind4_2 uwind5-uwind5_2 uwind6-uwind6_2 uwind7-uwind7_2 uwind8-uwind8_2 uwind9-uwind9_2 uwind10-uwind10_2 uwind11-uwind11_2];

resid=(tw_lhs-tw_rhs);

%-------------------------
%PLOTS

%Plot temperature cross-section
figure(1);
plot(xrange,temp_all);
xlabel('Latitude (deg. N)');
ylabel('Temperature (deg. C)');
title('Temperature cross-section on climatological 500 mb surface for January, 150 deg. E, 25-50 deg. N');

%Plot temperature gradient at each latitude
figure(2);
plot(xrange,temp_grad);
xlabel('Latitude (deg. N)');
ylabel('Temperature gradient (deg. C / meter)');
title('Temperature gradient cross-section on climatological 500 mb surface for January, 150 deg. E');

%Plot temperature gradient side of thermal wind
figure(3);
plot(xrange,tw_rhs);
xlabel('Latitude (deg. N)');
ylabel('Temperature gradient side of TW (m/s/mb)');
title('Temp grad side of TW for climatological 500 mb surface, January, 150 deg. E');

%Plot vertical wind shear side of thermal wind
figure(4);
plot(xrange,tw_lhs);
xlabel('Latitude (deg. N)');
ylabel('Vertical wind shear side of TW (m/s/mb)');
title('du/dp for climatological 500 mb surface, January, 150 deg. E');

%Plot residual
figure(5);
plot(xrange,resid);
xlabel('Latitude (deg. N)');
ylabel('Residual (m/s/mb)');
title('Residual of (vertical wind shear side - temperature gradient side) for climatological 500 mb surface, January, 150 deg. E');