% 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=640.*1000.;
drfine = 2000.;
pr = nrfine/drfine;

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

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

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

%x = [-1199.0*1000:2.0*1000:1199.0*1000];y = x;

load stretch_grid.dat
xx = stretch_grid;	%meters
yy = xx;
[x,y]=meshgrid(xx,yy);
%[x,y]=meshgrid(x,y);

rad = (x.^2+y.^2).^0.5;

%x = x*1000.;
%y = y*1000.;

% Get vorticity from subroutine of your choice:
for z=1:1    
vmax(z) = 0;
rmw(z) = 0;
vortfine=gaussianvort(rad,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

%whos vortfine
%fid = fopen('vorticity_sfc-nm02.dat','w')
%fwrite(fid,vortfine,'float')
%return;


%Vorticity Gradient
%[fx,fy] = gradient(vortfine,drfine);
[fx,fy] = gradient(vortfine,xx,yy);

output = fy';

fid = fopen('forcing_sfc-nm02.dat','w')
fwrite(fid,output,'float')

output = fx';

fid = fopen('forcing_sfc-nm02.dat','a')
fwrite(fid,output,'float')

return;
%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.dat','w');
%fwrite(fid,vt,'float');




%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')
