% Given the vorticity from a function, this subroutine
% computes the velocity profile

% Problem is that we need to integrate the vorticity
% which is defined on the stretched grid

% for this faster version, we do it once for a single radial
% grid that is higher resolution everywhere than the stretched grid.


%load new_vortex_parm.txt
%vmax = new_vortex_parm(:,1);        %m/s
%rmw = new_vortex_parm(:,2).*1000.;  %meters

load new_ht_grid.dat
z = new_ht_grid(:,1);	%meters

nrfine=300.*1000.;
drfine = 1000.;
pr = nrfine/drfine + 1;

vfine = zeros(pr,1);
vt = zeros(pr,71);

%drfine = r(nr)/(nrfine-1);

rfine = [0:drfine:nrfine];
[zbar,rbar]=meshgrid(z,rfine);

x = [-299:2.0:299];y = x;
[x,y]=meshgrid(x,y);
%x = x*1000.;
%y = y*1000.;

% Get vorticity from subroutine of your choice:
for z=1:71    
vmax(z) = 0.;
rmw(z) = 0.;
vortfine=gaussianvort(rfine,vmax(z),rmw(z));
%vortfine=vort3region(rfine);

%vortfine(1:30) = vortfine(30);


rvortfine=rfine.*vortfine;
vfine=drfine*cumtrapz(rvortfine);

vfine(2:pr)=vfine(2:pr)./rfine(2:pr);
%vfine(1)=0;

vt(:,z) = vfine';
vort(:,z) = vortfine';
end


%Parameters for tropical storm vortex (NM 02)
ztop = 14*1000;
alpha = 2.0;
c = 1.7;
gamma2 = 0;
delta2 = 0;
kappa = 0;
beta = 0;


% This formula is from Nolan and Montgomery 2002
vt = vt.*exp(-(c*zbar/ztop).^alpha);
%vt = vt.*exp(-(c*zbar/ztop).^alpha).* ...
%                (1 - gamma2*exp(-(rbar./(delta2 + kappa*zbar)).^beta) );


fid = fopen('ts-vortex-nm02-matlab-NEW.dat','w');
fwrite(fid,vt,'float');

return;




%Generate 2-D cylindrical grid
vt2 = zeros(pr,360,71);
for i=1:360
vt2(:,i,:) = vt;
end

%fid = fopen('stable_vortex_tan_fine.dat','w')
%fwrite(fid,vt,'float')

[r,t] = meshgrid(1.0:1.0:540.0,pi/180:pi/180:2*pi);
[nx,ny]=pol2cart(t,r);
fx = nx';fy = ny';


for k=1:71
z = vt2(:,:,k);
zi = griddata(fx,fy,z,x,y);
%nz(:,:,k) = zi;

%Converting to cartesian wind components
theta = atan2(y,x);
u = -1.*zi.*sin(theta);
u(isnan(u)) = 0;
v = zi.*cos(theta);
v(isnan(v)) = 0;

nu(:,:,k) = u';
nv(:,:,k) = v';
fprintf(1,'completed level:%i\n',k);
end

fid = fopen('forcing_uv71.dat','w')
fwrite(fid,nu,'float')
fid = fopen('forcing_uv71.dat','a')
fwrite(fid,nv,'float')

return;


%Vorticity
[curlz] = curl(x*1000.,y*1000.,u,v);

%Vorticity Gradient
[fx,fy] = gradient(curlz,drfine);

output = fy';

fid = fopen('forcing_uv71','w')
fwrite(fid,output,'double')

output = fx';

fid = fopen('forcing_uv71.dat','a')
fwrite(fid,output,'double')


%u = poicalc(fy,2,2,541,541);
%u = poicalc(single(fy),2,2,541,541);
return;

%Divergence
[div] = divergence(x*1000.,y*1000.,u,v);


%let u = if (y[g=nwind] ge 0) then (-1.*nwind*sin(the)) else (nwind*sin(the))
%let v = if (x[g=nwind] ge 0) then (nwind*cos(the)) else (-1.*nwind*cos(the))

%save('test.dat','vt(:,1)','-ASCII')
%fid = fopen('test_fine.dat','w')
%fwrite(fid,u','float')
%fid = fopen('test_fine.dat','a')
%fwrite(fid,v','float')
