% moisture advection at single pressure level from NCEP reanalysis data 
% 9-apr-08 AWR

% dependencies:
% NetCDF toolbox
% Dan Vimont's matlab toolbox:
% http://www.aos.wisc.edu/~dvimont/matlab.tar.gz
% mapping toolbox: http://www.eos.ubc.ca/~rich/map.html

clear
% global definitions for Dan Vimont's programs used to compute gradients
global RADUS RADIAN DEGREE
  if isempty(RADUS), RADUS=6375000; end
  if isempty(RADIAN), RADIAN=pi/180; end
  if isempty(DEGREE), DEGREE=180/pi; end;
  
%% get reanalysis daily u, v, q in Indian subdomain for an example year
%  SOURCES .NOAA .NCEP-NCAR .CDAS-1 .DAILY .Intrinsic .PressureLevel .u
%   P (850) VALUES
%   X (50E) (100E) RANGEEDGES
%   T (1 May 2001) (30 Sep 2001) RANGEEDGES
%   Y (30N) (0S) RANGEEDGES

% u-compt (m/s)
%loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP-NCAR/.CDAS-1/.DAILY/.Intrinsic/.PressureLevel/.u/P/%281000%29%28300%29RANGEEDGES/X/%2850E%29%28100E%29RANGEEDGES/T/%281%20May%202001%29%2830%20Sep%202001%29RANGEEDGES/Y/%2830N%29%280S%29RANGEEDGES/dods')

ncload('/Users/andy/data/India/Floods/data/u850_MaySep2001.nc');
% v-compt (m/s)
ncload('/Users/andy/data/India/Floods/data/v850_MaySep2001.nc');
% specific humidity (mass of water vapour per unit volume of air) (kg/kg?)
ncload('/Users/andy/data/India/Floods/data/qa850_MaySep2001.nc');
[nday nlev nlat nlon]=size(u);

%% moisture fluxes 
lev=1;
day=32;

ud=squeeze(u(day,lev,:,:));
vd=squeeze(v(day,lev,:,:));
qd=squeeze(qa(day,lev,:,:));

% moisture fluxes
uq=ud.*qd; vq=vd.*qd;

% moisture advection -(u.dq/dx + v.dq/dy)
[dqdx dqdy]=sph_grad(qd,X,Y);
udqdx=-ud.*dqdx; vdqdy=-vd.*dqdy;
adv=udqdx+vdqdy;

% moisture flux convergence
%[qdiv lat2 lon2]=sph_div1(uq,vq,X,Y);

%% plot maps to check
% winds
figure(1); clf
Xmat=repmat(X',length(Y),1);
Ymat=repmat(Y,1,length(X));
m_proj('Equidistant Cylindrical','longitudes',[50 100],'latitudes',[0 30]);
m_coast('patch',[.7 .7 .7],'edgecolor','none');
m_grid;
hold on
m_quiver(Xmat,Ymat,ud,vd,1,'color','black');
hold on
[cs,h]=m_contour(Xmat,Ymat,qd);
% clabel(cs,h,'manual');
title('uv & q fields 1 June 2001 (850hPa)')

%% fluxes
figure(2); clf
Xmat=repmat(X',length(Y),1);
Ymat=repmat(Y,1,length(X));
m_proj('Equidistant Cylindrical','longitudes',[50 100],'latitudes',[0 30]);
m_coast('patch',[.7 .7 .7],'edgecolor','none');
m_grid;
hold on
m_quiver(Xmat,Ymat,uq,vq,1,'color','black');
hold on
[cs,h]=m_contour(Xmat,Ymat,qd);
% clabel(cs,h,'manual');
title('moisture fluxes & q fields 1 June 2001 (850hPa)')

%% advection
figure(3); clf
Xmat=repmat(X',length(Y),1);
Ymat=repmat(Y,1,length(X));
m_proj('Equidistant Cylindrical','longitudes',[50 100],'latitudes',[0 30]);
m_coast('patch',[.7 .7 .7],'edgecolor','none');
m_grid;
hold on
m_quiver(Xmat,Ymat,udqdx,vdqdy,1,'color','black');
hold on
[cs,h]=m_contour(Xmat,Ymat,adv);
% clabel(cs,h,'manual');
title('moisture advection fields 1 June 2001 (850hPa)')

%print -f1 -dpdf name.pdf
