clear all;close all;clc

bryfile = 'FLC_v3_bry_SPINUP.nc';
grdfile = 'FLC_rossby_grd4_rx14spong3_MASK2.nc';
inifile='FLC_rossby_v3_ini.nc';
ncload(bryfile);

addpath /WSmounts/phoddat/share/domingues/MYROMS/matlab_NEWROMS/utility

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ncgrd = netcdf(grdfile);
h2 = ncgrd{'h'}(:,:);
lon = ncgrd{'lon_rho'}(:,:);
lat = ncgrd{'lat_rho'}(:,:);

lonv = ncgrd{'lon_v'}(:,:);
latv = ncgrd{'lon_v'}(:,:);

mask_rho = ncgrd{'mask_rho'}(:,:);
close(ncgrd);

N=28;
Vtransform=2;
Vstretching=4;
theta_s=4;
theta_b=0;
Tcline=10;

[zt]=set_depth(Vtransform, Vstretching, theta_s, theta_b, Tcline, N, 1, h2);
[ztV]=set_depth(Vtransform, Vstretching, theta_s, theta_b, Tcline, N, 3, h2);

lon_W = lon(:,1);
lat_W = lat(:,1);
mask_W = ~mask_rho(:,1)';


%  plot(lon_W,lat_W)
%  return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Defining inflow from FLC consitent with WOA13... adding barotropic component.

ncini = netcdf(inifile);
LAT_LIM = [22.7 24.8];

T_W = squeeze(ncini{'temp'}(1,:,:,1));
T_W2=T_W;
%  T_W(:,mask_W)=nan;
S_W = squeeze(ncini{'salt'}(1,:,:,1));
S_W2=S_W;
%  S_W(:,mask_W)=nan;


%  T_W = inpaint_nans(T_W,5);
%  S_W = inpaint_nans(S_W,5);

U_W = squeeze(ncini{'u'}(1,:,:,1));
V_W = squeeze(ncini{'v'}(1,:,:,1));

u.west = U_W';
v.west = V_W';

%  [ubar,vbar]=uv_barotropic(u,v,zt,[1 0 0 0]);

Z_W = squeeze(zt(:,1,:))';
lon_W = lon(:,1);
lat_W = lat(:,1);

K_FLC = find(lat_W >=LAT_LIM(1) & lat_W<=LAT_LIM(2));

T_FLC = T_W(:,K_FLC);
S_FLC = S_W(:,K_FLC);
Z_FLC = Z_W(:,K_FLC);
lat_FLC = lat_W(K_FLC);
lon_FLC = lon_W(K_FLC);

[lat2,aux] = meshgrid(lat_FLC,1:N);

figure
contourf(lat2,Z_FLC,T_FLC,0:30), shading flat, colorbar
xlim(LAT_LIM)
caxis([0 30])

figure
contourf(lat2,Z_FLC,S_FLC,30:.025:37), shading flat, colorbar
xlim(LAT_LIM)
caxis([34 38])

%  --------- Interpolate to Z 

Z_ref = 0:5:1450;

TZ_FLC = nan(length(Z_ref),length(lat_FLC));
SZ_FLC = TZ_FLC;

ZMAX = nanmax(-Z_FLC(:))

for i = 1:length(lat_FLC)

	sig = -Z_FLC(:,i);
	t_aux = T_FLC(:,i);
	s_aux = S_FLC(:,i);

	TZ_FLC(:,i) = interp1(sig,t_aux,Z_ref);
	SZ_FLC(:,i) = interp1(sig,s_aux,Z_ref);

end

MASK = isnan(TZ_FLC);

figure
pcolor(lat_FLC,-Z_ref,SZ_FLC), shading flat; colorbar
xlim(LAT_LIM)


U_FLC = get_vgeoSec(lon_FLC',lat_FLC',Z_ref,SZ_FLC,TZ_FLC,1449);
Knan=find(U_FLC>0);
U_FLC(Knan) = 0;

%  U_FLC = -U_FLC*3.0;
%  U_FLC = -U_FLC + 0.15;
U_FLC = -U_FLC;

U_FLC(MASK) = 0;

TSv_FLC = integrate_transpSec(lon_FLC',lat_FLC',Z_ref,U_FLC)

figure
pcolor(lat_FLC,-Z_ref,U_FLC), shading flat; colorbar
xlim(LAT_LIM)
polarcmap([-.7 .7],1)

UBAR_FLC = nanmean(U_FLC);

%  whos UBAR_FLC U_FLC lat_FLC
%  plot(lat_FLC,UBAR_FLC);
%  return
%  --------------------------------------------------------------------------
% interpolating to W boundary

for i = 1:length(lat_FLC)

	sig = -Z_FLC(:,i);
	u_aux = U_FLC(:,i);

	Usig_FLC(:,i) = interp1(Z_ref,u_aux,sig);

end

KU = find(isnan(Usig_FLC));
Usig_FLC(KU) = 0;

figure
contourf(lat2,Z_FLC,Usig_FLC,30), shading flat, colorbar
xlim(LAT_LIM)
polarcmap([-1.4 1.4],1)

figure
imagesc(U_W);colorbar
polarcmap([-1.4 1.4],1)

V_W(:,K_FLC) = Usig_FLC;
U_W(:,K_FLC) = Usig_FLC;

UBAR_W = zeros(1,length(lat_W));

UBAR_W(K_FLC) = UBAR_FLC;


figure
imagesc(V_W);colorbar
polarcmap([-1.4 1.4],1)
close all
keep S_W T_W U_W V_W clmfile T_USE W_WID Usig_FLC K_FLC T_FLC S_FLC UBAR_FLC grdfile ztV h2 v_north lonv latv UBAR_W
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  SETTING UP OUTFLOW ON N BRY

ncTSCLM = netcdf('FLC_rossby_v3_avg_tsCLM.nc');
v_north = squeeze(ncTSCLM{'v'}(1,:,209,:));
z_north = squeeze(ztV(209,:,:))';
close(ncTSCLM);

lon_north = lonv(209,:);
lat_north = latv(209,:);
lon_north2 = meshgrid(lonv(209,:),1:28);

K = find(v_north>100);
v_north(K) = 0;
Knan = isnan(v_north);
v_north(Knan) = 0;

FAC_v_north = abs(v_north)./nanmax(abs(v_north(:))); 

figure
contourf(lon_north2,z_north,FAC_v_north,-1:0.05:1), shading flat
caxis([-1 1])
colorbar

FAC_v_north = FAC_v_north*1.4;
v_north = v_north.*FAC_v_north;

figure
contourf(lon_north2,z_north,v_north,-1:0.05:1), shading flat
caxis([-.5 .5])
colorbar


z_ref = 0:5100;

for j=1:length(lon_north)
j
	v_aux = v_north(:,j);
	z_aux = z_north(:,j);

	Vz_north(:,j) = interp1(-z_aux,v_aux,z_ref);
end

Knan = find(isnan(Vz_north));
Vz_north(Knan) = 0;

figure
imagesc(Vz_north);
caxis([-.5 .5])
colorbar

TSv = integrate_transpSec(lon_north,lat_north,z_ref,Vz_north)
VBAR_NORTH = nanmean(Vz_north);

%whos lon_north v_north;
V_N = v_north;
imagesc(V_N);

close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  close all
nc = netcdf(bryfile,'write');

bry_time = nc{'bry_time'}(:)
%  bry_time(end)
bry_time = round(linspace(0,94608000,length(bry_time))) %3 years

u_west = nc{'u_west'}(:,:,:);
v_west = nc{'v_west'}(:,:,:);

ubar_west = nc{'ubar_west'}(:,:);
%  whos ubar_west UBAR_W;

v_north = nc{'v_north'}(:,:,:);
u_north = nc{'u_north'}(:,:,:);

vbar_north = nc{'vbar_north'}(:,:);

U_Naux = squeeze(u_north(1,:,:));
imagesc(U_Naux);

temp_west = nc{'temp_west'}(:,:,:);
salt_west = nc{'salt_west'}(:,:,:);

V_aux = squeeze(v_west(1,:,:));

for t=1:length(bry_time)

%  	u_west(t,:,:) = zeros(size(U_W));
%  	v_west(t,:,:) = V_W;

	if(t<13)
		UBAR_FAC = -0.02 + (t-1)*0.015;
	else
		UBAR_FAC = 0.16; 
	end

	[t UBAR_FAC]

	u_west(t,:,:) = U_W + UBAR_FAC;
	v_west(t,:,:) = zeros(size(V_aux));
	ubar_west(t,:) = UBAR_W + UBAR_FAC;

	v_north(t,:,:) = V_N;
	u_north(t,:,:) = zeros(size(U_Naux));

        vbar_north(t,:) = VBAR_NORTH;

	temp_west(t,:,:) = T_W2;
	salt_west(t,:,:) = S_W2;

end

NCWRITE=1;
if NCWRITE
	nc{'bry_time'}(:,:,:) = bry_time;

	nc{'u_west'}(:,:,:) = u_west;
	nc{'v_west'}(:,:,:) = v_west;

	nc{'ubar_west'}(:,:) = ubar_west;

	nc{'u_north'}(:,:,:) = u_north;
	nc{'v_north'}(:,:,:) = v_north;

        nc{'vbar_north'}(:,:) = vbar_north;

	nc{'temp_west'}(:,:,:) = temp_west;
	nc{'salt_west'}(:,:,:) = salt_west;

	result = sync(nc);


end
close(nc);
