% -------------------------------------------------- %
% Write All Input Files for MITgcm Offline Simulation
% Written By: Taylor Shropshire
% Date: 2/16/17
% -------------------------------------------------- %
clear;clc;close all;

year=2005;

% Paths 
% GoM Domain
path_main='/Net/gleam/taylor/GoM_50.1/';						   % all field dir 
path_in1=[path_main,num2str(year),'/netCDF/'];						   % year dynamical fields 
path_out1='/Net/gleam/taylor/postdoc/offline_fields/input_off/';			   % ouput
grid_file='/home/taylor/my_scripts/projects/calc_adv/GoM_domain/grid_GOMl0.04_03i.a';	   % hycom grid file 
bathy_file='/home/taylor/my_scripts/projects/calc_adv/GoM_domain/depth_GOMl0.04_03i.a';   % hycom bathy file
river_file='/Net/gleam/taylor/GoM_50.1/rivers/rivers_03i.a';				   % hycom river file 

% Functions 
addpath /home/taylor/my_scripts/projects/MITgcm						% lldistkm
addpath /nexsan/people/taylor/projects/hycom_bb86/BB86_PACKAGE/MATLAB/UTILITIES		% subvar2

% Load variables
load('/nexsan/people/taylor/projects/NPZ/data/matlab/eastcoastmask.mat','cc')

% Original GoM Domain (Parameters in HYCOM native .a archive files to netCDF conversion script )
idm=541;	% x dimension of full field
jdm=385;	% y dimension of fill field 

% Script flags
write_intial=0;			% write inital input files (bathy etc.)
write_out_flag=0;

% NetCDF file names
nc_fnames_uv_tran=dir([path_in1,'uv_tran/*.nc']);		% u and v transports
nc_fnames_ts=dir([path_in1,'ts_fields/*.nc']);			% t and s fields

% ---------------------- Inital File Processing --------------------- % 
display('Reading Inital File and Preprocessing'); display(''); 

tt=1;
nc_file1=[path_in1,'ts_fields/',nc_fnames_ts(tt).name];	ncid1=netcdf.open(nc_file1,'NOWRITE');
nc_file2=[path_in1,'uv_tran/',nc_fnames_uv_tran(tt).name]; ncid2=netcdf.open(nc_file2,'NOWRITE');
varid=netcdf.inqVarID(ncid1,'Latitude'); lat=netcdf.getVar(ncid1,varid); lat=double(lat);
varid=netcdf.inqVarID(ncid1,'Longitude'); lon=netcdf.getVar(ncid1,varid); lon=double(lon);	
varid=netcdf.inqVarID(ncid1,'Depth'); depth_uv=netcdf.getVar(ncid1,varid); depth_uv=double(depth_uv);
varid=netcdf.inqVarID(ncid1,'pot_temp'); tmp=netcdf.getVar(ncid1,varid); tmp=double(tmp);
varid=netcdf.inqVarID(ncid1,'salinity'); sal=netcdf.getVar(ncid1,varid); sal=double(tmp);
varid=netcdf.inqVarID(ncid2,'ut'); u_tran = netcdf.getVar(ncid2,varid);	u_tran=double(u_tran);
varid=netcdf.inqVarID(ncid2,'vt'); v_tran = netcdf.getVar(ncid2,varid);	v_tran=double(v_tran);
netcdf.close(ncid1);netcdf.close(ncid2);

% Get w interfaces
depth_w=0;
for zz=1:length(depth_uv)
depth_w(zz+1,1)=depth_w(zz)+(depth_uv(zz)-depth_w(zz))*2;
end 
dz=diff(depth_w);

% Load grid descriptors and bathy 
plon=sub_var2(grid_file,[jdm,idm],19,1);        % Lon of the T point 
plat=sub_var2(grid_file,[jdm,idm],19,2);        % Lat of the T point
gdx=sub_var2(grid_file,[jdm,idm],19,16);        % dx of t cell (centered on v variable)
gdy=sub_var2(grid_file,[jdm,idm],19,15);        % dy of t cell (centered on u variable)
bathy_hycom=sub_var2(bathy_file,[jdm,idm],1,1);
% Note: bathy_hycom is not used to generate the MITgcm bathy file 
% intead it is generated from the velocity field to ensure that masks agree (see below).

% Make grid from netCDF coordaintes
[lon,lat]=meshgrid(lon,lat);

% Rotate arrays
u_tran2=[]; v_tran2=[]; sal2=[]; tmp2=[];
for zz=1:size(tmp,3);
u_tran2(:,:,zz)=u_tran(:,:,zz)'; 	v_tran2(:,:,zz)=v_tran(:,:,zz)';
tmp2(:,:,zz)=tmp(:,:,zz)';		sal2(:,:,zz)=sal(:,:,zz)';
end
u_tran=u_tran2; clear u_tran2;		v_tran=v_tran2; clear v_tran2;
tmp=tmp2; clear tmp2;			sal=sal2; clear sal2;
% Note: variables read in from netCDF need to be rotated but variables read in by sub_var2 dont

% Adjust dimension for c-grid. 
plon=plon(1:end-1,1:end-1);	plat=plat(1:end-1,1:end-1);
tmp=tmp(1:end-1,1:end-1,:);	sal=sal(1:end-1,1:end-1,:);
bathy_hycom=bathy_hycom(1:end-1,1:end-1);
xdim_ori=size(bathy_hycom,2);    ydim_ori=size(bathy_hycom,1);
u_tran=u_tran(1:end-1,:,:); 	gdy=gdy(1:end-1,:);
v_tran=v_tran(:,1:end-1,:);	gdx=gdx(:,1:end-1);

% Set Fill Values
fill_value=max(reshape(tmp,[1,numel(tmp)]));
u_tran(find(u_tran==fill_value | u_tran==0))=NaN; 	% must set 0 to NaN for mask to correspond 
v_tran(find(v_tran==fill_value | v_tran==0))=NaN;	% must set 0 to NaN for mask to correspond 
tmp(find(tmp==fill_value))=NaN;
sal(find(sal==fill_value))=NaN;
bathy_hycom(find(bathy_hycom==fill_value))=NaN;

% **** Optional subset (indicies for t cell) **** % 
%x1=1; 	y1=1;
%x2=540; y2=384;

x1=1; x2=size(bathy_hycom,2);
y1=1; y2=size(bathy_hycom,1);

plon=plon(y1:y2,x1:x2);         	plat=plat(y1:y2,x1:x2);
tmp=tmp(y1:y2,x1:x2,:);                 sal=sal(y1:y2,x1:x2,:);
u_tran=u_tran(y1:y2,x1:x2+1,:);		gdy=gdy(y1:y2,x1:x2+1); 	
v_tran=v_tran(y1:y2+1,x1:x2,:);		gdx=gdx(y1:y2+1,x1:x2); 	
bathy_hycom=bathy_hycom(y1:y2,x1:x2);
xdim=size(tmp,2);	ydim=size(tmp,1);
xc=plon;		yc=plat;
% ********************************************* % 

% --------- % Adjust bathy for z-Level Full Cells ------ % 
bathy_comp=[];			% used to compare to bathy field generated from dynamical fields (below)
for yy=1:size(bathy_hycom,1)
for xx=1:size(bathy_hycom,2)

if ~isnan(bathy_hycom(yy,xx))
ind=find(depth_w>bathy_hycom(yy,xx));

if isempty(ind)
bathy_comp(yy,xx)=depth_w(end);
else
bathy_comp(yy,xx)=depth_w(ind(1));
end

else 
bathy_comp(yy,xx)=bathy_hycom(yy,xx);
end

end; end
% -------------------------------------------% 

% ----------- Create Masks ----------------- % 
w_mask1=ones(size(bathy_comp,1),size(bathy_comp,2),length(depth_w));
for xx=1:size(bathy_comp,2)
for yy=1:size(bathy_comp,1)
if isnan(bathy_comp(yy,xx))==1
w_mask1(yy,xx,:)=NaN;
else
ind=find(depth_w>=bathy_comp(yy,xx));        w_mask1(yy,xx,ind)=NaN;
end
end
end

display('Creating Masks'); display('')
u_mask=u_tran;		u_mask(find(~isnan(u_mask)))=1;
v_mask=v_tran;		v_mask(find(~isnan(v_mask)))=1;
t_mask=tmp;		t_mask(find(~isnan(t_mask)))=1;
s_mask=sal;		s_mask(find(~isnan(s_mask)))=1;

w_mask2=ones(size(bathy_comp,1),size(bathy_comp,2),length(depth_w));
for zz=1:length(depth_uv)
for xx=1:size(bathy_comp,2)
for yy=1:size(bathy_comp,1)
u1m=u_mask(yy,xx,zz); u2m=u_mask(yy,xx+1,zz);
v1m=v_mask(yy,xx,zz); v2m=v_mask(yy+1,xx,zz);
temp=nanmean([u1m u2m v1m v2m]); 	% all cell faces must be NaN if the cell is land

if isnan(temp)==1
w_mask2(yy,xx,zz)=NaN;
else
w_mask2(yy,xx,zz)=1;
end

end; end; end
w_mask2(:,:,end)=NaN;

% Check mask agreement.
ind1=find(isnan(w_mask1));
ind2=find(isnan(w_mask2));	% 162 additional cells than ind1 (but includes all of ind1 - necessary)
ind3=setdiff(ind1,ind2);
if ~isempty(ind3)
display('Mask Error'); pause 
end 

% Set mask 
w_mask=w_mask2; clear w_mask1 w_mask2; 

% Must create new bathy from u-vel and v-vel to get mask for T and S (order is important)
bathy=[];
for xx=1:size(w_mask,2)
for yy=1:size(w_mask,1)
ind=find(isnan(w_mask(yy,xx,:)));
bathy(yy,xx)=depth_w(ind(1));
end
end
bathy(find(bathy==0))=NaN;

% Create T and S mask
ts_mask=ones(size(bathy,1),size(bathy,2),length(depth_uv));
for xx=1:size(bathy,2)
for yy=1:size(bathy,1)

if isnan(bathy(yy,xx))==1
ts_mask(yy,xx,:)=NaN;
else
ind=find(depth_uv>=bathy(yy,xx)); 	ts_mask(yy,xx,ind)=NaN;
end

end
end

display('')
display('Finished Inital File Readings');
display('*****************************');
display('Begin Main File Read')
display('')

% ####################################################################
% ######## Read NC Files for Flow Fields and Write Input Files #######
% ####################################################################
u_lim=[0,0]; v_lim=[0,0]; w_lim=[0,0]; tmp_lim=[0,0]; sal_lim=[0,0];
u_mask_cells=[]; v_mask_cells=[]; t_mask_cells =[]; s_mask_cells=[]; 
file_count=1;

% *********** Main Loop ************** %
record_start=1;
record_end=length(nc_fnames_uv_tran);
for tt=record_start:record_end
tic
display(nc_fnames_uv_tran(tt).name)
uv_file=[path_in1,'uv_tran/',nc_fnames_uv_tran(tt).name]; 	ncid1=netcdf.open(uv_file,'NOWRITE');
ts_file=[path_in1,'ts_fields/',nc_fnames_ts(tt).name];		ncid2=netcdf.open(ts_file,'NOWRITE');
ncid1=netcdf.open(uv_file,'NOWRITE');	ncid2=netcdf.open(ts_file,'NOWRITE');
varid=netcdf.inqVarID(ncid1,'ut'); u_tran = netcdf.getVar(ncid1,varid);		u_tran=double(u_tran);
varid=netcdf.inqVarID(ncid1,'vt'); v_tran = netcdf.getVar(ncid1,varid);		v_tran=double(v_tran);
varid=netcdf.inqVarID(ncid2,'pot_temp'); tmp = netcdf.getVar(ncid2,varid);	tmp=double(tmp);
varid=netcdf.inqVarID(ncid2,'salinity'); sal = netcdf.getVar(ncid2,varid);	sal=double(sal);
netcdf.close(ncid1); netcdf.close(ncid2);

% Rotate Array
u_tran2=[]; v_tran2=[]; tmp2=[]; sal2=[];
for zz=1:size(u_tran,3);
u_tran2(:,:,zz)=u_tran(:,:,zz)';	v_tran2(:,:,zz)=v_tran(:,:,zz)';
tmp2(:,:,zz)=tmp(:,:,zz)';		sal2(:,:,zz)=sal(:,:,zz)';
end
u_tran=u_tran2; clear u_tran2;		v_tran=v_tran2; clear v_tran2;
tmp=tmp2; clear tmp2;			sal=sal2; clear sal2;

% Get variables on c-grid 
u_tran=u_tran(1:end-1,:,:); v_tran=v_tran(:,1:end-1,:);
tmp=tmp(1:end-1,1:end-1,:); sal=sal(1:end-1,1:end-1,:);

% Remove Fill Values 
u_tran(find(u_tran==fill_value | u_tran==0))=0;		% set to zero here instead of NaN (done above)
v_tran(find(v_tran==fill_value | v_tran==0))=0;		% set to zero here instead of NaN (done above)
tmp(find(tmp==fill_value))=NaN;	
sal(find(sal==fill_value))=NaN;

% --- Record land cells ---- %
u_mask_cells(1,tt)=length(find(u_tran==0));
v_mask_cells(1,tt)=length(find(v_tran==0));
t_mask_cells(1,tt)=length(find(isnan(tmp)));
s_mask_cells(1,tt)=length(find(isnan(sal)));
% --------------------------- % 

% Optional subset - full field if x1=1, x2=xdim
u_tran=u_tran(y1:y2,x1:x2+1,:);		v_tran=v_tran(y1:y2+1,x1:x2,:);
tmp=tmp(y1:y2,x1:x2,:);			sal=sal(y1:y2,x1:x2,:);

% ---- New ---- % 
% ts_mask is applied to tmp and sal below
u_tran=u_tran.*u_mask; u_tran(find(isnan(u_tran)))=0;
v_tran=v_tran.*v_mask; v_tran(find(isnan(v_tran)))=0;
% -------------- % 

% ******************* Fix TS Field ***************************** % 

% ------------- Neariest (non zero) neighbor interpolation ------------ % 
[x_field,y_field]=meshgrid([1:size(tmp,2)],[1:size(tmp,1)]);
num_elements=size(tmp,1)*size(tmp,2);

tmp(find(isnan(tmp)))=9999;      tmp=tmp.*ts_mask; 	tmp2=tmp;
sal(find(isnan(sal)))=9999;      sal=sal.*ts_mask;	sal2=sal;	

% Interpolation for each level
for zz=1:size(tmp,3)
t_surf=tmp(:,:,zz);                             s_surf=sal(:,:,zz);
[ind_bady,ind_badx]=find(s_surf==9999);

% Interpolation for each bad point 
if ~isempty(ind_bady) & zz<26                  % don't do interp for deep layers (takes too much time)
for ii=1:length(ind_bady);
xc_point=xc(ind_bady(ii),ind_badx(ii));        % lon value
yc_point=yc(ind_bady(ii),ind_badx(ii));        % lat value 

% Extract a small box around the bad point
SD=1;
for nn=1:num_elements
ind1=find(x_field<=ind_badx(ii)+SD & x_field>=ind_badx(ii)-SD & y_field<=ind_bady(ii)+SD & y_field>=ind_bady(ii)-SD);
xc_sample=xc(ind1);             yc_sample=yc(ind1);
t_sample=t_surf(ind1);          s_sample=s_surf(ind1);

t_sample(find(t_sample==9999))=NaN;
ind2=find(isnan(t_sample));
xc_sample(ind2)=[];             yc_sample(ind2)=[];
t_sample(ind2)=[];              s_sample(ind2)=[];

if isempty(t_sample)    % get bigger if nothing is found
SD=SD+SD;
else
break
end

end % nn

% Calculate distance to possible NN points and replace bad point
dist=[];
for jj=1:length(t_sample)
dist(jj)=lldistkm([yc_point,xc_point],[yc_sample(jj),xc_sample(jj)]);
end
ind3=find(dist==min(dist));
t_replace=mean(t_sample(ind3));         % ind3 can be more than 1 
s_replace=mean(s_sample(ind3));         % ind3 can be more than 1 

tmp2(ind_bady(ii),ind_badx(ii),zz)=t_replace;
sal2(ind_bady(ii),ind_badx(ii),zz)=s_replace;
end % ii
% Note: could save ind3 indicies and not have to calculate again for each field

else    % lower level correction (prety sure this takes hours using NN interp so I just use avg)
ind=find(t_surf~=9999 & ~isnan(t_surf));
t_replace=mean(t_surf(ind));	s_replace=mean(s_surf(ind));
for ii=1:length(ind_bady)
tmp2(ind_bady(ii),ind_badx(ii),zz)=t_replace;
sal2(ind_bady(ii),ind_badx(ii),zz)=s_replace;
end % ii

end % if ~isempty(ind_bady) & zz<26
end % zz

% Update TS field and check that all bad cells were fixed
tmp=tmp2;clear tmp2; 	sal=sal2;clear sal2;
ind=find(tmp==9999 | sal==9999);
ind1=find(isnan(ts_mask));		ind2=find(isnan(tmp));
if ~isempty(ind) | length(ind1)~=length(ind2)
display('TS Adjustment Error'); pause
end
% *********************** END TS mask Adjustment ********************* % 

% ************************** Compute W Field ********************** % 

% Get velocities from transport
u_vel=[];
for zz=1:length(depth_w)-1
for xx=1:size(u_tran,2)
for yy=1:size(u_tran,1)
ay=gdy(yy,xx)*dz(zz);
u_vel(yy,xx,zz)=u_tran(yy,xx,zz)/ay;
end
end
end

v_vel=[];
for zz=1:length(depth_w)-1
for xx=1:size(v_tran,2)
for yy=1:size(v_tran,1)
ax=gdx(yy,xx)*dz(zz);
v_vel(yy,xx,zz)=v_tran(yy,xx,zz)/ax;
end
end
end

% Adjust flow fields so that domain discretization can be done with multiple processors
x_new1=1;	x_new2=536;
y_new1=2;	y_new2=381;
% Note: first row (1,:) has bogus velocities always start with at least 2

% Trim for parallel processing 
u_vel=u_vel(y_new1:y_new2,x_new1:x_new2+1,:);
v_vel=v_vel(y_new1:y_new2+1,x_new1:x_new2,:);
sal=sal(y_new1:y_new2,x_new1:x_new2,:);
tmp=tmp(y_new1:y_new2,x_new1:x_new2,:);

if tt==record_start
gdx2=gdx(y_new1:y_new2+1,x_new1:x_new2,:);
gdy2=gdy(y_new1:y_new2,x_new1:x_new2+1,:);
u_mean=u_vel*0; v_mean=v_vel*0;
w_mean=zeros(size(tmp,1),size(tmp,2),size(tmp,3)+1);
sal_mean=sal*0; tmp_mean=tmp*0;
xc2=xc(y_new1:y_new2,x_new1:x_new2); yc2=yc(y_new1:y_new2,x_new1:x_new2);
whos u_vel v_vel tmp sal
end 

% Calculate w-vel
w_vel=zeros(size(tmp,1),size(tmp,2),size(tmp,3)+1);

for zz=length(depth_w)-1:-1:1
for xx=1:size(tmp,2)
for yy=1:size(tmp,1)

dx=gdx2(yy,xx);		% changes with lat
dx2=gdx2(yy+1,xx);	% changes with lat
dy=gdy2(yy,xx);         % dy and dy2 should be the same since dy does not change with longitude
dy2=gdy2(yy,xx+1);	% dy and dy2 should be the same since dy does not change with longitude

w_vel(yy,xx,zz)=2*dz(zz)/(dy*(dx+dx2))*...
((u_vel(yy,xx,zz)*dy-u_vel(yy,xx+1,zz)*dy2)+...
(v_vel(yy,xx,zz)*dx-v_vel(yy+1,xx,zz)*dx2))+...
w_vel(yy,xx,zz+1);

end
end
end

% Mask check
temp_mask1=w_mask(y_new1:y_new2,x_new1:x_new2,:);
temp_mask2=w_vel;					temp_mask2(find(temp_mask2==0))=NaN;
ind1=find(isnan(temp_mask1));                         	ind2=find(isnan(temp_mask2));
ind3=setdiff(ind1,ind2);                        	% finds indices that are in ind1 that are not in ind2 
if ~isempty(ind3)
display('Mask Error'); pause
end
% ************************** End Compute W-Field ******************** % 

% Find Maximum and Minimum Values
temp=reshape(u_vel,[1,numel(u_vel)]); 	temp(find(temp==0))=[];
u_max=max(temp);u_min=min(temp);
temp=reshape(v_vel,[1,numel(v_vel)]);	temp(find(temp==0))=[];
v_max=max(temp);v_min=min(temp);
temp=reshape(w_vel,[1,numel(w_vel)]);	temp(find(temp==0))=[];
w_max=max(temp);w_min=min(temp);
temp=reshape(tmp,[1,numel(tmp)]);		temp(find(temp==0))=[];
tmp_max=max(temp);tmp_min=min(temp);
temp=reshape(sal,[1,numel(sal)]);		temp(find(temp==0))=[];
sal_max=max(temp);sal_min=min(temp);

if tt==record_start
u_lim(1,1)=u_min; u_lim(1,2)=u_max;
v_lim(1,1)=v_min; v_lim(1,2)=v_max;
w_lim(1,1)=w_min; w_lim(1,2)=w_max;
tmp_lim(1,1)=tmp_min; tmp_lim(1,2)=tmp_max;
sal_lim(1,1)=sal_min; sal_lim(1,2)=sal_max;
else
u_lim(1,1)=min([u_lim(1,1),u_min]);u_lim(1,2)=max([u_lim(1,2),u_max]);
v_lim(1,1)=min([v_lim(1,1),v_min]);v_lim(1,2)=max([v_lim(1,2),v_max]);
w_lim(1,1)=min([w_lim(1,1),w_min]);w_lim(1,2)=max([w_lim(1,2),w_max]);
tmp_lim(1,1)=min([tmp_lim(1,1),tmp_min]);tmp_lim(1,2)=max([tmp_lim(1,2),tmp_max]);
sal_lim(1,1)=min([sal_lim(1,1),sal_min]);sal_lim(1,2)=max([sal_lim(1,2),sal_max]);
end

display(['u(Min,Max) = ','(',num2str(u_lim(1,1)),',',num2str(u_lim(1,2)),')']);
display(['v(Min,Max) = ','(',num2str(v_lim(1,1)),',',num2str(v_lim(1,2)),')']);
display(['w(Min,Max) = ','(',num2str(w_lim(1,1)),',',num2str(w_lim(1,2)),')']);
display(['tmp(Min,Max) = ','(',num2str(tmp_lim(1,1)),',',num2str(tmp_lim(1,2)),')']);
display(['sal(Min,Max) = ','(',num2str(sal_lim(1,1)),',',num2str(sal_lim(1,2)),')']);

% Sum to compute average 
u_mean=u_mean+u_vel;		v_mean=v_mean+v_vel;
tmp_mean=tmp_mean+tmp;		sal_mean=sal_mean+sal;
w_mean=w_mean+w_vel;

% **************** Write Flow Field Files ******************
if write_out_flag==1
display('Writting Output Files')

% Set all arrays to the same dimensions for writting
u_vel=u_vel(:,1:end-1,:);		v_vel=v_vel(1:end-1,:,:);
tmp(find(isnan(tmp)))=0;		sal(find(isnan(sal)))=0;
w_vel=w_vel(:,:,1:end-1);

% Check for NaN
if ~isempty(find(isnan(tmp))) | ~isempty(find(isnan(sal))) | ~isempty(find(isnan(u_vel))) | ~isempty(find(isnan(v_vel))) ...
| ~isempty(find(isnan(w_vel)))
display('NaN Error before writing');pause
end

% Rotate for Writting
u_vel2=[]; v_vel2=[]; w_vel2=[]; tmp2=[]; sal2=[];
for zz=1:size(u_vel,3)
u_vel2(:,:,zz)=u_vel(:,:,zz)'; 	v_vel2(:,:,zz)=v_vel(:,:,zz)';
tmp2(:,:,zz)=tmp(:,:,zz)';	sal2(:,:,zz)=sal(:,:,zz)';
w_vel2(:,:,zz)=w_vel(:,:,zz)';
end
u_vel=u_vel2;clear u_vel2;	v_vel=v_vel2;clear v_vel2;
tmp=tmp2;clear tmp2;		sal=sal2;clear sal2;
w_vel=w_vel2;clear w_vel2;

% Write Daily Files
file_name1=[path_out1,'uvel/u.',num2str(tt,'%4.4i'),'.data'];
file_name2=[path_out1,'vvel/v.',num2str(tt,'%4.4i'),'.data'];
file_name3=[path_out1,'wvel/w.',num2str(tt,'%4.4i'),'.data'];
file_name4=[path_out1,'tmp/t.',num2str(tt,'%4.4i'),'.data'];
file_name5=[path_out1,'sal/s.',num2str(tt,'%4.4i'),'.data'];
fid1=fopen(file_name1,'w','b'); fwrite(fid1,u_vel,'float32'); fclose(fid1);
fid2=fopen(file_name2,'w','b'); fwrite(fid2,v_vel,'float32'); fclose(fid2);
fid3=fopen(file_name3,'w','b'); fwrite(fid3,w_vel,'float32'); fclose(fid3);
fid4=fopen(file_name4,'w','b'); fwrite(fid4,tmp,'float32'); fclose(fid4);
fid5=fopen(file_name5,'w','b'); fwrite(fid5,sal,'float32'); fclose(fid5);
else
display('Not Writting Output Files')
end 
% ***** End Writting Flow Field Files ****** % 

display('Done one file');pause

file_count=file_count+1;
display(''); toc
end % tt

display('Finished Writting Flow Field Files')
pause
% #############################################################################
% ############### End Reading NC Files for Flow Fields ########################
% #############################################################################

% Generate average flow field figures 
u_mean=u_mean/tt;	u_mean(find(u_mean==0))=NaN;
v_mean=v_mean/tt;	v_mean(find(v_mean==0))=NaN;
w_mean=w_mean/tt;	w_mean(find(w_mean==0))=NaN;
tmp_mean=tmp_mean/tt;	tmp_mean(find(tmp_mean==0))=NaN;
sal_mean=sal_mean/tt;	sal_mean(find(sal_mean==0))=NaN;

figure(1); pcolor(u_mean(:,:,1));shading flat;colorbar; colormap(jet); caxis([-1 1])
figure(2); pcolor(v_mean(:,:,1));shading flat;colorbar; colormap(jet); caxis([-1 1])
figure(3); pcolor(w_mean(:,:,1));shading flat;colorbar; colormap(jet); caxis([-1e-6 1e-6])
figure(4); pcolor(tmp_mean(:,:,1));shading flat;colorbar; colormap(jet); caxis([20 30])
figure(5); pcolor(sal_mean(:,:,1));shading flat;colorbar; colormap(jet); caxis([33 37])

% ###### Write Inital Condition, Boundary Condition, Bathy, Relax, Grid Fiels ########## %
if write_inital==1

% Trim mask
ts_mask=ts_mask(y_new1:y_new2,x_new1:x_new2,:);

% ************ Write Bathy ************ %
bathy=bathy(y_new1:y_new2,x_new1:x_new2);
bathy=bathy*-1;					% Negative for Mitgcm 	
bathy=bathy'; 					% Rotate back for writting as is done with advection fields
bathy(find(isnan(bathy)))=0;
file_name4=[path_out2,'depth_gom.bin'];		fid1=fopen(file_name4,'w','b');
fwrite(fid1,bathy,'float32');			fclose(fid1);
% ************************************* %

% ***************** River *************** % 
dx_rA=sub_var2(grid_file,[jdm,idm],19,16);
dy_rA=sub_var2(grid_file,[jdm,idm],19,15);
dx_rA=dx_rA(y1_ori:y2_ori,x1_ori:x2_ori);
dy_rA=dy_rA(y1_ori:y2_ori,x1_ori:x2_ori);
dx_rA=dx_rA(:,1:end-1);				dy_rA=dy_rA(1:end-1,:);
dx_rA=dx_rA(y1:y2+1,x1:x2);			dy_rA=dy_rA(y1:y2,x1:x2+1);

% USGS river data daily climatology
file='/home/taylor/my_scripts/projects/USGS_river/daily_climatology_nitrate_simplified.xlsx';
data_all=importdata(file);		data=data_all.data;
NO_r=data(:,end);			NO_r=NO_r*1/1000*1000/62.0049*1000;

rA=[];
for xx=1:size(dx_rA,2)
for yy=1:size(dx_rA,1)-1
rA(yy,xx)=(dx_rA(yy+1,xx)+dx_rA(yy,xx))/2*dy_rA(yy,xx);
end
end
surf_vol=rA.*dz(1);

riv_dis=[]; riv_nutr=[];
for ii=1:12
riv=sub_var2(river_file,[jdm,idm],12,ii);
riv=riv(y1_ori:y2_ori,x1_ori:x2_ori);	riv=riv(1:end-1,1:end-1);
riv=riv(y1:y2,x1:x2);			riv_dis(:,:,ii)=riv.*rA;
end

% Define time stamp for monthly discharge
month_ind=[1,32,60,91,121,152,182,213,244,274,305,335,366];  % First day of the month for each month
time=[1:365]; month_mid_ind=[];
for mm=1:length(month_ind)-1;
month_mid_ind(mm)=(month_ind(mm)+month_ind(mm+1))/2;
end

d1_field=[]; d1_field=riv_dis(:,:,end);d1_field(:,:,2)=riv_dis(:,:,1);
t1=[month_mid_ind(end),365+month_mid_ind(1)];
t2=[month_mid_ind(end):0.5:365+month_mid_ind(1)];

% Interpolate to first and last day
field=abs(riv_dis);field=sum(field,3); d3=[]; 
for xx=1:size(riv_dis,2)
for yy=1:size(riv_dis,1)
if field(yy,xx)>0
d1=squeeze(d1_field(yy,xx,:))';	d2=interp1(t1,d1,t2);  ind=find(t2==366);
d3(yy,xx)=d2(ind);
else
d3(yy,xx)=0;
end; end; end
% why not just d3(yy,xx)=interp1(t1,d1,366);

% Add first and last fields so that daily interpolation can be done.
riv_dis2=[];
riv_dis2(:,:,2:size(riv_dis,3)+1)=riv_dis;
riv_dis2(:,:,1)=d3; 				riv_dis2(:,:,end+1)=d3;
riv_dis=riv_dis2;clear riv_dis2;
month_new=[1,month_mid_ind,366];

NO_riv=NO_r; 	NO_riv(60)=[];  % February 29th 

% Interpolate discharge to all days in year
riv_dis2=zeros(size(riv_dis,1),size(riv_dis,2),length(time));
riv_nutr=riv_dis2;
for xx=1:size(riv_dis,2)
for yy=1:size(riv_dis,1)
if field(yy,xx)>0
temp=squeeze(riv_dis(yy,xx,:));
riv_dis2(yy,xx,:)=interp1(month_new,temp,time);		riv_nutr(yy,xx,:)=NO_riv;
end
end
end
riv_dis=riv_dis2; clear riv_dis2;

% Check
field2=abs(riv_dis);field2=sum(field2,3);
ind1=find(field~=0);	ind2=find(field2~=0);
if length(ind1)~=length(ind2)
display('Error');pause;
end

% --- New --- % 
for tt=1:size(riv_nutr,3);
surf=riv_nutr(:,:,tt); surf(ind_land)=0; riv_nutr(:,:,tt)=surf;
surf=riv_dis(:,:,tt); surf(ind_land)=0; riv_dis(:,:,tt)=surf;
end
% ------------- %

% Trim for parallel processing
riv_nutr=riv_nutr(y_new1:y_new2,x_new1:x_new2,:);
riv_dis=riv_dis(y_new1:y_new2,x_new1:x_new2,:);
ind1=find(riv_nutr~=0);

riv_nutr2=[];riv_dis2=[];
for tt=1:size(riv_nutr,3)
riv_nutr2(:,:,tt)=riv_nutr(:,:,tt)'.*ts_mask(:,:,1)';
riv_dis2(:,:,tt)=riv_dis(:,:,tt)'.*ts_mask(:,:,1)';
end
riv_nutr=riv_nutr2;clear riv_nutr2; 	riv_dis=riv_dis2;clear riv_dis2;
riv_nutr(find(isnan(riv_nutr)))=0;	riv_dis(find(isnan(riv_dis)))=0;
ind2=find(riv_nutr~=0);

% Check
if length(ind1)~=length(ind2)
display('River Field Error'); pause
end

% Write 
file_name5=[path_out2,'riv_nutr_daily_080820_OB.bin'];
file_name6=[path_out2,'riv_dis_daily_080820_OB.bin'];
fid1=fopen(file_name5,'w','b');	fwrite(fid1,riv_nutr,'float32'); fclose(fid1);
fid2=fopen(file_name6,'w','b');	fwrite(fid2,riv_dis,'float32'); fclose(fid2);
% *********************************************************** % 

% *********************** Light Data **************** % 
% -------- CFSR (1993-2010) -------- % 
path_in2='/Net/gleam/taylor/postdoc/light_forcing/dom_080820/day_avg/';
nc_fnames_light=dir([path_in2,'*TaqaQrQp_GOMSAB.nc']);

% Get lat lon for CFSR fields 
ii=1;
file=([path_in2,nc_fnames_light(ii).name]); ncid=netcdf.open(file,'NOWRITE');
varid=netcdf.inqVarID(ncid,'Latitude');lat_light=netcdf.getVar(ncid,varid); lat_light=double(lat_light);
varid=netcdf.inqVarID(ncid,'Longitude');lon_light=netcdf.getVar(ncid,varid); lon_light=double(lon_light);
netcdf.close(ncid);

% Make grid 
[lon_lg,lat_lg]=meshgrid(lon_light,lat_light);

% Process All Years
light_count=1; 		light_all=[];
%for ii=1:length(nc_fnames_light)
%ii=13;					% 2005
ii=1;
nc_fnames_light(ii).name
file=([path_in2,nc_fnames_light(ii).name]);	ncid=netcdf.open(file,'NOWRITE');
varid=netcdf.inqVarID(ncid,'shwflx');rad=netcdf.getVar(ncid,varid); rad=double(rad);
netcdf.close(ncid);

albed=0.09;		% Value used by the Navy - Alex Bozec
rad=rad*(1-albed);
% Note: Don't use if working with net shortwave radiation

% Interpolate to model grid 
for tt=1:size(rad,3)
surf=rad(:,:,tt)';			surf2=interp2(lon_lg,lat_lg,surf,xc,yc);
light_all(:,:,light_count)=surf2;	light_count=light_count+1;
end

% Trim for parallel processing
light_all=light_all(y_new1:y_new2,x_new1:x_new2,:);

% Apply mask and process for writting
light_all2=[];
for tt=1:size(light_all,3)
light_all2(:,:,tt)=light_all(:,:,tt)'.*ts_mask(:,:,1)';
end
light_all=light_all2;clear light_all2;

% Write Light Forcing File
light_all(find(isnan(light_all)))=0;
file_name5=[path_out2,'light_data_080820_OB.bin'];
fid1=fopen(file_name5,'w','b');	fwrite(fid1,light_all,'float32'); fclose(fid1);
% *************************************************** % 

% *********************** Inital Condtions **************** % 
load('/nexsan/people/taylor/projects/mit_gcm/simulation_projects/1_grid_point_GoM/GoM1_120718_expt001/details/input/bgc_prof_120618.mat')

% Took off one level in 1D model runs
NO_prof(end+1)=NO_prof(end);		NH_prof(end+1)=NH_prof(end);
SI_prof(end+1)=SI_prof(end);		DON_prof(end+1)=DON_prof(end);
PON_prof(end+1)=PON_prof(end);		OP_prof(end+1)=OP_prof(end);
sp_prof(end+1)=sp_prof(end);		lp_prof(end+1)=lp_prof(end);
sz_prof(end+1)=sz_prof(end);		lz_prof(end+1)=lz_prof(end);
pz_prof(end+1)=pz_prof(end);

% Rotate
NO=[];NH=[];SI=[];DON=[];PON=[];OP=[];sp=[];lp=[];sz=[];lz=[];pz=[];
for zz=1:size(ts_mask,3)
temp=ts_mask(:,:,zz)';
NO(:,:,zz)=temp*NO_prof(zz); 	NH(:,:,zz)=temp*NH_prof(zz);
SI(:,:,zz)=temp*SI_prof(zz);	DON(:,:,zz)=temp*DON_prof(zz);
PON(:,:,zz)=temp*PON_prof(zz);	OP(:,:,zz)=temp*OP_prof(zz);
sp(:,:,zz)=temp*sp_prof(zz);	lp(:,:,zz)=temp*lp_prof(zz);
sz(:,:,zz)=temp*sz_prof(zz);	lz(:,:,zz)=temp*lz_prof(zz);
pz(:,:,zz)=temp*pz_prof(zz);
end

% Remove NaN;
ind=find(isnan(NO));
NO(ind)=0;	NH(ind)=0;	SI(ind)=0;
DON(ind)=0;	PON(ind)=0;	OP(ind)=0;
sp(ind)=0;	lp(ind)=0;	sz(ind)=0;
lz(ind)=0;	pz(ind)=0;

% Write files
file_name1=[path_out2,'NO_field.bin'];		file_name2=[path_out2,'NH_field.bin'];
file_name3=[path_out2,'SI_field.bin'];		file_name4=[path_out2,'DON_field.bin'];
file_name5=[path_out2,'PON_field.bin'];		file_name6=[path_out2,'OP_field.bin'];
file_name7=[path_out2,'sp_field.bin'];		file_name8=[path_out2,'lp_field.bin'];
file_name9=[path_out2,'sz_field.bin'];		file_name10=[path_out2,'lz_field.bin'];
file_name11=[path_out2,'pz_field.bin'];
fid1=fopen(file_name1,'w','b'); fwrite(fid1,NO,'float32'); fclose(fid1);
fid2=fopen(file_name2,'w','b'); fwrite(fid2,NH,'float32'); fclose(fid2);
fid3=fopen(file_name3,'w','b'); fwrite(fid3,SI,'float32'); fclose(fid3);
fid4=fopen(file_name4,'w','b'); fwrite(fid4,DON,'float32'); fclose(fid4);
fid5=fopen(file_name5,'w','b'); fwrite(fid5,PON,'float32'); fclose(fid5);
fid6=fopen(file_name6,'w','b'); fwrite(fid6,OP,'float32'); fclose(fid6);
fid7=fopen(file_name7,'w','b'); fwrite(fid7,sp,'float32'); fclose(fid7);
fid8=fopen(file_name8,'w','b'); fwrite(fid8,lp,'float32'); fclose(fid8);
fid9=fopen(file_name9,'w','b'); fwrite(fid9,sz,'float32'); fclose(fid9);
fid10=fopen(file_name10,'w','b'); fwrite(fid10,lz,'float32'); fclose(fid10);
fid11=fopen(file_name11,'w','b'); fwrite(fid11,pz,'float32'); fclose(fid11);
% ******************************************* % 

% ************ Open Boundary Condtions ************************ %
% Northern Boundary
OBN_NO=squeeze(NO(:,end,:)); 		OBN_NH=squeeze(NH(:,end,:)); 
OBN_SI=squeeze(SI(:,end,:)); 		OBN_DON=squeeze(DON(:,end,:));
OBN_PON=squeeze(PON(:,end,:));		OBN_OP=squeeze(OP(:,end,:)); 
OBN_sp=squeeze(sp(:,end,:)); 		OBN_lp=squeeze(lp(:,end,:)); 
OBN_sz=squeeze(sz(:,end,:)); 		OBN_lz=squeeze(lz(:,end,:)); 
OBN_pz=squeeze(pz(:,end,:)); 
file_name1=[path_out2,'OBN_NO.bin'];		file_name2=[path_out2,'OBN_NH.bin'];
file_name3=[path_out2,'OBN_SI.bin'];		file_name4=[path_out2,'OBN_DON.bin'];
file_name5=[path_out2,'OBN_PON.bin'];		file_name6=[path_out2,'OBN_OP.bin'];
file_name7=[path_out2,'OBN_sp.bin'];		file_name8=[path_out2,'OBN_lp.bin'];
file_name9=[path_out2,'OBN_sz.bin'];		file_name10=[path_out2,'OBN_lz.bin'];
file_name11=[path_out2,'OBN_pz.bin'];
fid1=fopen(file_name1,'w','b'); fwrite(fid1,OBN_NO,'float32'); fclose(fid1);
fid2=fopen(file_name2,'w','b'); fwrite(fid2,OBN_NH,'float32'); fclose(fid2);
fid3=fopen(file_name3,'w','b'); fwrite(fid3,OBN_SI,'float32'); fclose(fid3);
fid4=fopen(file_name4,'w','b'); fwrite(fid4,OBN_DON,'float32'); fclose(fid4);
fid5=fopen(file_name5,'w','b'); fwrite(fid5,OBN_PON,'float32'); fclose(fid5);
fid6=fopen(file_name6,'w','b'); fwrite(fid6,OBN_OP,'float32'); fclose(fid6);
fid7=fopen(file_name7,'w','b'); fwrite(fid7,OBN_sp,'float32'); fclose(fid7);
fid8=fopen(file_name8,'w','b'); fwrite(fid8,OBN_lp,'float32'); fclose(fid8);
fid9=fopen(file_name9,'w','b'); fwrite(fid9,OBN_sz,'float32'); fclose(fid9);
fid10=fopen(file_name10,'w','b'); fwrite(fid10,OBN_lz,'float32'); fclose(fid10);
fid11=fopen(file_name11,'w','b'); fwrite(fid11,OBN_pz,'float32'); fclose(fid11);

% Eastern Boundary
OBE_NO=squeeze(NO(end,:,:));   		OBE_NH=squeeze(NH(end,:,:));   
OBE_SI=squeeze(SI(end,:,:));   		OBE_DON=squeeze(DON(end,:,:)); 
OBE_PON=squeeze(PON(end,:,:)); 		OBE_OP=squeeze(OP(end,:,:));   
OBE_sp=squeeze(sp(end,:,:));   		OBE_lp=squeeze(lp(end,:,:));   
OBE_sz=squeeze(sz(end,:,:));   		OBE_lz=squeeze(lz(end,:,:));   
OBE_pz=squeeze(pz(end,:,:));   
file_name1=[path_out2,'OBE_NO.bin'];		file_name2=[path_out2,'OBE_NH.bin'];
file_name3=[path_out2,'OBE_SI.bin'];		file_name4=[path_out2,'OBE_DON.bin'];
file_name5=[path_out2,'OBE_PON.bin'];		file_name6=[path_out2,'OBE_OP.bin'];
file_name7=[path_out2,'OBE_sp.bin'];		file_name8=[path_out2,'OBE_lp.bin'];
file_name9=[path_out2,'OBE_sz.bin'];		file_name10=[path_out2,'OBE_lz.bin'];
file_name11=[path_out2,'OBE_pz.bin'];
fid1=fopen(file_name1,'w','b'); fwrite(fid1,OBE_NO,'float32'); fclose(fid1);
fid2=fopen(file_name2,'w','b'); fwrite(fid2,OBE_NH,'float32'); fclose(fid2);
fid3=fopen(file_name3,'w','b'); fwrite(fid3,OBE_SI,'float32'); fclose(fid3);
fid4=fopen(file_name4,'w','b'); fwrite(fid4,OBE_DON,'float32'); fclose(fid4);
fid5=fopen(file_name5,'w','b'); fwrite(fid5,OBE_PON,'float32'); fclose(fid5);
fid6=fopen(file_name6,'w','b'); fwrite(fid6,OBE_OP,'float32'); fclose(fid6);
fid7=fopen(file_name7,'w','b'); fwrite(fid7,OBE_sp,'float32'); fclose(fid7);
fid8=fopen(file_name8,'w','b'); fwrite(fid8,OBE_lp,'float32'); fclose(fid8);
fid9=fopen(file_name9,'w','b'); fwrite(fid9,OBE_sz,'float32'); fclose(fid9);
fid10=fopen(file_name10,'w','b'); fwrite(fid10,OBE_lz,'float32'); fclose(fid10);
fid11=fopen(file_name11,'w','b'); fwrite(fid11,OBE_pz,'float32'); fclose(fid11);

% Southern Boundary
OBS_NO=squeeze(NO(:,1,:));   		OBS_NH=squeeze(NH(:,1,:));   
OBS_SI=squeeze(SI(:,1,:));   		OBS_DON=squeeze(DON(:,1,:)); 
OBS_PON=squeeze(PON(:,1,:)); 		OBS_OP=squeeze(OP(:,1,:));   
OBS_sp=squeeze(sp(:,1,:));   		OBS_lp=squeeze(lp(:,1,:));   
OBS_sz=squeeze(sz(:,1,:));   		OBS_lz=squeeze(lz(:,1,:));   
OBS_pz=squeeze(pz(:,1,:));   
file_name1=[path_out2,'OBS_NO.bin'];		file_name2=[path_out2,'OBS_NH.bin'];
file_name3=[path_out2,'OBS_SI.bin'];		file_name4=[path_out2,'OBS_DON.bin'];
file_name5=[path_out2,'OBS_PON.bin'];		file_name6=[path_out2,'OBS_OP.bin'];	
file_name7=[path_out2,'OBS_sp.bin'];		file_name8=[path_out2,'OBS_lp.bin'];
file_name9=[path_out2,'OBS_sz.bin'];		file_name10=[path_out2,'OBS_lz.bin'];
file_name11=[path_out2,'OBS_pz.bin'];
fid1=fopen(file_name1,'w','b'); fwrite(fid1,OBS_NO,'float32'); fclose(fid1);
fid2=fopen(file_name2,'w','b'); fwrite(fid2,OBS_NH,'float32'); fclose(fid2);
fid3=fopen(file_name3,'w','b'); fwrite(fid3,OBS_SI,'float32'); fclose(fid3);
fid4=fopen(file_name4,'w','b'); fwrite(fid4,OBS_DON,'float32'); fclose(fid4);
fid5=fopen(file_name5,'w','b'); fwrite(fid5,OBS_PON,'float32'); fclose(fid5);
fid6=fopen(file_name6,'w','b'); fwrite(fid6,OBS_OP,'float32'); fclose(fid6);
fid7=fopen(file_name7,'w','b'); fwrite(fid7,OBS_sp,'float32'); fclose(fid7);
fid8=fopen(file_name8,'w','b'); fwrite(fid8,OBS_lp,'float32'); fclose(fid8);
fid9=fopen(file_name9,'w','b'); fwrite(fid9,OBS_sz,'float32'); fclose(fid9);
fid10=fopen(file_name10,'w','b'); fwrite(fid10,OBS_lz,'float32'); fclose(fid10);
fid11=fopen(file_name11,'w','b'); fwrite(fid11,OBS_pz,'float32'); fclose(fid11);
% ****************************************************** % 

% ***************** Grid Descriptors ********************* % 
% Grid descriptiors must bt +1 on both dim 
x2_ori=x2_ori+1;	y2_ori=y2_ori+1;
x2=x2+1;		y2+y2+1;
x_new2=x_new2+1;	y_new2=y_new2+1;
% can't run code above now that this is changed

% XC YC (tracer points) 
plon=sub_var2(grid_file,[jdm,idm],19,1);
plat=sub_var2(grid_file,[jdm,idm],19,2);
plon(find(plon>180))=plon(find(plon>180))-360;
plon=plon(y1_ori:y2_ori,x1_ori:x2_ori);		plat=plat(y1_ori:y2_ori,x1_ori:x2_ori);
plon=plon(1:end-1,1:end-1);			plat=plat(1:end-1,1:end-1);
plon=plon(y1:y2,x1:x2);				plat=plat(y1:y2,x1:x2);
XC=plon(y_new1:y_new2,x_new1:x_new2);		YC=plat(y_new1:y_new2,x_new1:x_new2);

% XG YG (vorticity point)
qlon=sub_var2(grid_file,[jdm,idm],19,3);
qlat=sub_var2(grid_file,[jdm,idm],19,4);
qlon(find(qlon>180))=qlon(find(qlon>180))-360;
qlon=qlon(y1_ori:y2_ori,x1_ori:x2_ori);     	qlat=qlat(y1_ori:y2_ori,x1_ori:x2_ori);
qlon=qlon(1:end-1,1:end-1);                     qlat=qlat(1:end-1,1:end-1);
qlon=qlon(y1:y2,x1:x2);				qlat=qlat(y1:y2,x1:x2);
XG=qlon(y_new1:y_new2,x_new1:x_new2);           YG=qlat(y_new1:y_new2,x_new1:x_new2);

% dxF and dyF (Dim of Tracer Cell at the Center)  
pcsx=sub_var2(grid_file,[jdm,idm],19,10);
pcsy=sub_var2(grid_file,[jdm,idm],19,11);
pcsx=pcsx(y1_ori:y2_ori,x1_ori:x2_ori);         pcsy=pcsy(y1_ori:y2_ori,x1_ori:x2_ori);
pcsx=pcsx(:,1:end-1);				pcsy=pcsy(1:end-1,:);
pcsx=pcsx(y1:y2+1,x1:x2);			pcsy=pcsy(y1:y2,x1:x2+1);
dxF=pcsx(y_new1:y_new2+1,x_new1:x_new2);	dyF=pcsy(y_new1:y_new2,x_new1:x_new2+1);

% rA (Dim of Tracer Cell for Area Calculation)
vcsx=sub_var2(grid_file,[jdm,idm],19,16);
ucsy=sub_var2(grid_file,[jdm,idm],19,15);
vcsx=vcsx(y1_ori:y2_ori,x1_ori:x2_ori);		ucsy=ucsy(y1_ori:y2_ori,x1_ori:x2_ori);
vcsx=vcsx(:,1:end-1);                           ucsy=ucsy(1:end-1,:);
vcsx=vcsx(y1:y2+1,x1:x2);                       ucsy=ucsy(y1:y2,x1:x2+1);
dx_rA=vcsx(y_new1:y_new2+1,x_new1:x_new2);      dy_rA=ucsy(y_new1:y_new2,x_new1:x_new2+1);

% dxV dyU (Dim of Q Cell at Center)
qcsx=sub_var2(grid_file,[jdm,idm],19,12);
qcsy=sub_var2(grid_file,[jdm,idm],19,13);
qcsx=qcsx(y1_ori:y2_ori,x1_ori:x2_ori);         qcsy=qcsy(y1_ori:y2_ori,x1_ori:x2_ori);
qcsx=qcsx(:,1:end-1);                           qcsy=qcsy(1:end-1,:);
qcsx=qcsx(y1:y2+1,x1:x2);                       qcsy=qcsy(y1:y2,x1:x2+1);
dxV=qcsx(y_new1:y_new2+1,x_new1:x_new2);        dyU=qcsy(y_new1:y_new2,x_new1:x_new2+1);

% rAz (Dim of Q Cell for Area Calculation)
ucsx=sub_var2(grid_file,[jdm,idm],19,14);
vcsy=sub_var2(grid_file,[jdm,idm],19,17);
ucsx=ucsx(y1_ori:y2_ori,x1_ori:x2_ori);         vcsy=vcsy(y1_ori:y2_ori,x1_ori:x2_ori);
ucsx=ucsx(:,1:end-1);                           vcsy=vcsy(1:end-1,:);
ucsx=ucsx(y1:y2+1,x1:x2);                       vcsy=vcsy(y1:y2,x1:x2+1);
dx_rAz=ucsx(y_new1:y_new2+1,x_new1:x_new2);     dy_rAz=vcsy(y_new1:y_new2,x_new1:x_new2+1);

% Area of T Cell
rA=[];rAz=[];rAw=[];rAs=[];
for xx=1:size(dx_rA,2)
for yy=1:size(dx_rA,1)-1
rA(yy,xx)=(dx_rA(yy+1,xx)+dx_rA(yy,xx))/2*dy_rA(yy,xx);         % Area of T Cell
rAz(yy,xx)=(dx_rAz(yy+1,xx)+dx_rAz(yy,xx))/2*dy_rAz(yy,xx);     % Area of Q Cell
rAw(yy,xx)=(dxV(yy+1,xx)+dxV(yy,xx))/2*dyF(yy,xx);              % Area of U Cell
rAs(yy,xx)=(dxF(yy+1,xx)+dxF(yy,xx))/2*dyU(yy,xx);              % Area of V Cell
end
end
% Note: Since dx changes with latitude for a single cell (Trapizoid shape) the area is x1+x2/2*y

% All grid descriptors must have the same array size in file
dxC=dx_rAz(1:end-1,:);		dyC=dy_rAz(:,1:end-1);
dxG=dx_rA(1:end-1,:);		dyG=dy_rA(:,1:end-1);
dxV=dxV(1:end-1,:);		dyF=dyF(:,1:end-1);
dxF=dxF(1:end-1,:);		dyU=dyU(:,1:end-1);

% Display dimensions
whos tmp XC YC dxF dyF rA XG YG dxV dyU rAz dxC dyC rAw rAs dxG dyG

grid=zeros(size(XC,2),size(XC,1),18);
grid(:,:,1)=XC';				grid(:,:,2)=YC';
grid(:,:,3)=dxF';				grid(:,:,4)=dyF';
grid(:,:,5)=rA';				grid(:,:,6)=XG';
grid(:,:,7)=YG';				grid(:,:,8)=dxV';
grid(:,:,9)=dyU';				grid(:,:,10)=rAz';
grid(:,:,11)=dxC';				grid(:,:,12)=dyC';
grid(:,:,13)=rAw';				grid(:,:,14)=rAs';
grid(:,:,15)=dxG';				grid(:,:,16)=dyG';
grid(:,:,17)=dyG'*0+1;    
grid(:,:,18)=dyG'*0;      
fid=fopen([path_out2,'grid.face001.bin'],'w','b'); fwrite(fid,grid,'double'); fclose(fid);

end 

