% ------------------------------------------------------------------------- %
%  Read MITgcm netCDF Output Fiels 
% Written by: Taylor Shropshire
% Date: July 2018
% Description: Reads in each individual day across all tiles 
% then organizes tiles into one 3D field. Diagnostics are computed 
% for each record and saved. Previous code would read all days and then 
% compute diagnostics. This was fine for test simulations with smaller amounts of data
% ------------------------------------------------------------------------- %
clear;clc;close all

% Paths 
path1='/Net/ftp_pub/taylor/GoM1x_expt001_111919/';
path2='/nexsan/people/taylor/projects/NPZ/data/matlab/eastcoastmask.mat';
path3='/nexsan/people/taylor/projects/mit_gcm/MITgcm2/MITgcm/verification/gom_sim/output_fig/';
path8='/home/taylor/my_scripts/projects/calc_adv/GoM_domain/depth_GOMl0.04_03i.a';

% Load
load(path2,'cc');	% east coast mask

% Get all tile directory names
dir_names=dir([path1,'/*tile*']); 

% ################################################### % 
% ############### Model Parameters ################## % 
% ################################################### % 

% ************ Non Bio Parameters ************* % 
t_step=60*60;                   % time step (s)
dt_b=t_step/86400;              % time step used in code (1/day)
par_frac=0.43;                  % fraction of photosynthetically active radiation
sink=15;                        % sinking of PON and OP
vert_diff=1.75e-4;               % Vertical diffusivity (not used in the code)
min_val=1.0e-6;                 % min value 

sed_frac=0.0;                   % fraction of bottom cell sinking that goes towards sedimentation
sig_mid=0.2;                    % sigmodial function for sedimentaiton
sig_slope=40;                   % slope of sigmodial functon for sedimentation
sed_max=0.0;                    % maximum sedimentation 
n2c_conv=106/16*1/1000*12.0107*1000; % mmol N to mg C
% ************ END Non Bio Parameters *********** % 

% ********************** Bio Parameters ************************* %
% Light Atteuation
par_frac=0.43;                          % fraction of photosynthetically active radiation
ext_w=0.03;                             % attenuation due to water - [0.04]
ext_sp = 0.03;                          % atteunation due to small phytoplankton - [0.04]
ext_lp =0.03;                           % atteunation due to large phytoplankton - [0.04] 

% Temperature Dependence
tdep1=0.0693; scale_tmp_phyto1=1;        % phytoplankton growth
tdep2=0.0693; scale_tmp_phyto2=1;        % phytoplankton mortality
tdep3=0.0519; scale_tmp_phyto3=1;        % phytolpankton respiration
tdep4=0.0693; scale_tmp_zoo1=1;          % zooplankton grazing (1:0.0693, 0.14:0.1386, 0.016:0.2079)
tdep5=0.0693; scale_tmp_zoo2=1;          % zooplankton mortality (1:0.0693, 0.14:0.1386)
tdep6=0.0693; scale_tmp_decomp=1;        % decomposition (1:0.0693, 0.14:0.1386)

% Small Phytoplankton 
vmax_sp = 0.4;                          % max growth rate at 0 deg C (1/d) - [0.4]
k_NO_sp = 0.5;                          % nitrate half saturation constant - [1.0]
k_NH_sp = 0.1;                          % ammonium half saturation constant - [0.1]
alpha_sp = 0.1;                         % PI curve parameter alpha - [0.01]
beta_sp = 4.5e-4;                       % PI curve parmater beta - [4.5e-4], 1.4e-3
ref_resp_sp = 0.03;                     % respiration at 0 deg C (1/d) - [0.03]
ref_mort_sp = 0.002;                    % mortality at at 0 deg C (1/d) - [0.0585], Linear Mort - ~[0.002]
ext_excr_sp = 0.135;                    % extracellular excretion - [0.135]
tc_v_sp = tdep1;                        % tmp coefficent for growth - [0.0693]
tc_r_sp = tdep3;                        % tmp coefficent for respiration - [0.0519]
tc_m_sp = tdep2;                        % tmp coeffiect for mortality - [0.0693]
inh_NH_NO_sp = 0;                       % nitrate uptake inhibiton, NOTE: parameter currently not being used - [1.4] 
chl2c_sp_const=1/50;                 % chlorophyll : carbon ratio
sp_mort_type = 1;                       % mortality type

% Large Phytoplankton 
vmax_lp = 0.8;                          % max growth rate at 0 deg C (1/d) - [0.8]
k_NO_lp = 3.0;                          % nitrate half saturation constant - [3.0]
k_NH_lp = 0.3;                          % ammonium half saturation constant - [0.3]
k_SI_lp = 6.0;                          % silica half saturation constant - [6.0] 
alpha_lp = 0.1;                         % PI curve parameter alpha - [0.01], 1.4e-3
beta_lp= 4.5e-4;                        % PI curve parmater beta - [4.5e-4]
ref_resp_lp = 0.03;                     % respiration at 0 deg C (1/d) - [0.03]
ref_mort_lp = 0.001;                    % mortality at at 0 deg C (1/d) - [0.029], Linear Mort - ~[0.001]
ext_excr_lp = 0.135;                    % extracellular excretion - [0.0135]
tc_v_lp = tdep1;                        % tmp coefficent for growth - [0.0693]
tc_r_lp = tdep3;                        % tmp coefficent for respiration - [0.0519]
tc_m_lp = tdep2;                        % tmp coeffiect for mortality - [0.0693]
inh_NH_NO_lp = 0;                       % nitrate uptake inhibiton, NOTE: parameter currently not being used - [1.4]
chl2c_lp_const=0.03;                 % chlorophyll : carbon ratio
lp_mort_type = 1/25;                       % mortality type 

% Small Zooplankton 
gmax_sz_sp = 0.6;                       % max grazing rate on small phytoplankton at 0 deg C (1/d) - [0.04]
gmax2_sz_sp = 0.0;                      % minumum grazing under no light. (~0.1)
rad_graz_sz_sp = 0.0;                   % half saturation constant for light limited grazing (~10.0)     
tc_g_sz_sp = tdep4;                     % tmp coefficent for grazing - [0.0693]
iv_sz_sp = 1.4;                         % ivlev constant - [1.4]
thresh_sz_sp = 0.04;                    % feeding threshold on small phytoplankton - [0.043]
tc_m_sz = tdep5;                        % tmp coefficent for mortality - [0.0693]
ae_sz = 0.7;                            % assimilation efficency - [0.7]
gge_sz = 0.3;                           % gross growth efficency - [0.3]
ref_mort_sz = 0.022;                     % mortality at 0 deg C (1/d) - [0.0585], Linear Mort - ~[0.002]
sz_mort_type = 1;                       % mortality type 

% Large Zooplankton 
gmax_lz_sp = 0.0;                       % max grazing rate on small phytoplankton at 0 deg C (1/d) - [0.1]
gmax_lz_lp = 0.3;                       % max grazing rate on large phytoplankton at 0 deg C (1/d) - [0.4]
gmax_lz_sz = 0.3;                       % max grazing rate on small zooplankton at 0 deg C (1/d) - [0.4]
tc_g_lz_sp = tdep4;                     % tmp coefficent for grazing -[0.0693]
tc_g_lz_lp = tdep4;                     % tmp coefficent for grazing -[0.0693]
tc_g_lz_sz = tdep4;                     % tmp coefficent for grazing -[0.0693]
iv_lz_sp = 1.4;                         % ivlev constant - [1.4]
iv_lz_lp = 1.4;                         % ivlev constant - [1.4]
iv_lz_sz = 1.4;                         % ivlev constant - [1.4]
thresh_lz_sp = 100.0;                   % feeding threshold on small phytoplankton - [0.04]
thresh_lz_lp = 0.04;                    % feeding threshold on large phytoplankton - [0.04]
thresh_lz_sz = 0.04;                    % feeding threshold on small zooplankton  - [0.04]
ref_mort_lz = 0.022;                     % mortality at at 0 deg C (1/d) - [0.0585], Linear Mort - ~[0.002]
tc_m_lz = tdep5;                        % tmp coefficent for mortality - [0.0693]
ae_lz = 0.7;                            % assimilation efficency - [0.7]
gge_lz = 0.3;                           % gross growth efficency - [0.3]
lz_mort_type = 1;                       % mortality type

% Predatory Zooplankton 
gmax_pz_lp = 0.1;                       % max grazing rate on large phytoplankton at 0 deg C (1/d) - [0.2]
gmax_pz_sz = 0.1;                       % max grazing rate on small zooplankton at 0 deg C (1/d) - [0.2]
gmax_pz_lz = 0.3;                       % max grazing rate on large zooplankton at 0 deg C (1/d) - [0.2]
tc_g_pz_lp = tdep4;                     % tmp coefficent for grazing -[0.0693]
tc_g_pz_sz = tdep4;                     % tmp coefficent for grazing -[0.0693]
tc_g_pz_lz = tdep4;                     % tmp coefficent for grazing -[0.0693]
iv_pz_lp = 1.4;                         % ivlev constant - [1.4]
iv_pz_sz = 1.4;                         % ivlev constant - [1.4]
iv_pz_lz = 1.4;                         % ivlev constant - [1.4]
thresh_pz_lp = 0.04;                    % feeding threshold on large phytoplankton - [0.04]
thresh_pz_sz = 0.04;                    % feeding threshold on small zooplankton - [0.04]
thresh_pz_lz = 0.04;                    % feeding threshold on large zooplankton - [0.04]
inh_szlz_lp = 4.605;                    % grazing on large phytoplankton inhibition by small and large zooplankton - [4.605]
inh_lz_sz = 3.01;                       % grazing on small zooplankton inhibition by large zooplankton - [3.01]
ref_mort_pz = 0.12;                     % mortality at at 0 deg C (1/d) - [0.0585], Linear Mort - ~[0.002]
tc_m_pz = tdep5;                        % tmp coefficent for mortality - [0.0693]
ae_pz = 0.7;                            % assimilation efficency - [0.7]
gge_pz = 0.3;                           % gross growth efficency - [0.3] 
pz_mort_type = 2;                       % mortality type 

% Nutrients (11)
ref_nitr = 0.003;                       % nitrification at 0 deg C (1/d) - [0.03]
ref_dec_PON_NH = 0.01;                  % decompositon from PON to NH at 0 deg C (1/d) - [0.1]
ref_dec_PON_DON = 0.05;                 % decompositon from PON to DON at 0 deg C (1/d) - [0.1]
ref_dec_DON_NH = 0.02;                  % decompositon from DON to NH at 0 deg C (1/d) - [0.2] 
ref_dec_OP_SI = 0.01;                   % decompositon from OP to SI at 0 deg C (1/d) - [0.1]
r_SI_N = 1.0;                           % ratio of silica to nitrogen (SI/N) - [2.0]
tc_nitr = tdep6;                        % tmp coefficent for nitrification - [0.0693] 
tc_dec_PON_NH = tdep6;                  % tmp coefficent for decomposition  - [0.0693]
tc_dec_PON_DON = tdep6;                 % tmp coefficent for decomposition - [0.0693]
tc_dec_DON_NH = tdep6;                  % tmp coefficent for decomposition - [0.0693]
tc_dec_OP_SI = tdep6;                   % tmp coefficent for decomposition - [0.0693]

min_val=1.0e-6;                         % min value (used in the Silica Limitation/Nitrogen Limitation Term for LP

% Chl2c Submodel (Li et al., 2010)
chl2c_sp_min = 0.0000;                  % minimum Chl:C ratio - [0.000]
chl2c_lp_min = 0.005;                   % maximum Chl:C ratio - [0.005] 
chl2c_sp_max = 0.015;                   % minimum Chl:C ratio - [0.03]
chl2c_lp_max = 0.03;                    % maximum Chl:C ratio - [0.061] 
alpha_chl_sp = 0.28;                    % chlorophyll specific initial slope of PE curve - [0.28 +- 0.11]
alpha_chl_lp = 0.28;                    % chlorophyll specific initial slope of PE curve - [0.28 +- 0.11]
% ############# END Bio Parameters ################### %

% ======== Get grid coordinates and spatial dimension of output files ======== % 
xc_all=[];yc_all=[];rA_all=[];bathy_all=[];cfrac_all=[];

% ------ Read each tile -------- %  
for dd=1:length(dir_names)
grid_file=dir([path1,dir_names(dd).name,'/*grid*']);
nc_file1=[path1,dir_names(dd).name,'/',grid_file(1).name];ncid1=netcdf.open(nc_file1,'NOWRITE');

if dd==1
varid=netcdf.inqVarID(ncid1,'RC');depth_uv=netcdf.getVar(ncid1,varid);depth_uv=depth_uv'*-1;
varid=netcdf.inqVarID(ncid1,'RF');depth_w=netcdf.getVar(ncid1,varid);depth_w=depth_w*-1;
varid=netcdf.inqVarID(ncid1,'drF');depth_dz=netcdf.getVar(ncid1,varid);depth_dz=depth_dz;
end

varid=netcdf.inqVarID(ncid1,'XC');xc=netcdf.getVar(ncid1,varid);
varid=netcdf.inqVarID(ncid1,'YC');yc=netcdf.getVar(ncid1,varid);
varid=netcdf.inqVarID(ncid1,'Depth');bathy=netcdf.getVar(ncid1,varid);
varid=netcdf.inqVarID(ncid1,'HFacC');cfrac=netcdf.getVar(ncid1,varid);
varid=netcdf.inqVarID(ncid1,'rA');rA=netcdf.getVar(ncid1,varid);
netcdf.close(ncid1);

% Rotate 
xc=xc';yc=yc';rA=rA';bathy=bathy';
cfrac2=[];
for zz=1:size(cfrac,3)
cfrac2(:,:,zz)=cfrac(:,:,zz)';
end
cfrac=cfrac2;clear cfrac2;

% Keep all tiles (organize below into a single field)
xc_all(:,:,dd)=xc; yc_all(:,:,dd)=yc;
bathy_all(:,:,dd)=bathy; cfrac_all(:,:,:,dd)=cfrac;
rA_all(:,:,dd)=rA;
end % dd
% ----------------------------- %

% Define domain lat lon limits
temp1=reshape(xc_all,[1,numel(xc_all)]);
temp2=reshape(yc_all,[1,numel(yc_all)]);
dom_x1=min(temp1); dom_x2=max(temp1);
dom_y1=min(temp2); dom_y2=max(temp2);

% Define all tile lat lon limits
xlim_all=[]; ylim_all=[];
xlim_all(:,1)=squeeze(xc_all(1,1,:));	% bottom left x coordinate for all tiles 
xlim_all(:,2)=squeeze(xc_all(1,end,:));	% bottom right x corrdinate for all tiles
ylim_all(:,1)=squeeze(yc_all(1,1,:));	% bottom left y coordinate for all tiles
ylim_all(:,2)=squeeze(yc_all(end,1,:));	% top left y coorodinate for all tiles

% Define unique tile lat lon limits 
xlim=[]; ylim=[];
xlim(:,1)=sort(unique(xlim_all(:,1)));xlim(:,2)=sort(unique(xlim_all(:,2)));
ylim(:,1)=sort(unique(ylim_all(:,1)));ylim(:,2)=sort(unique(ylim_all(:,2)));

% Dimensions
txd=size(xc,2); tyd=size(yc,1); tzd=length(depth_uv);		% define tile dimension
num_tx=size(xlim,1); num_ty=size(ylim,1);			% number of tiles in x and y

% ----- Organize into one array ------ %
xc_all2=[];yc_all2=[];rA_all2=[];bathy_all2=[];cfrac_all2=[];
y1=1; y2=tyd;
for yy=1:num_ty
x1=1; x2=txd;

for xx=1:num_tx
ind=find(xlim_all(:,1)==xlim(xx,1) & xlim_all(:,2)==xlim(xx,2)...
& ylim_all(:,1)==ylim(yy,1) & ylim_all(:,2)==ylim(yy,2));

xc_all2(y1:y2,x1:x2)=xc_all(:,:,ind);yc_all2(y1:y2,x1:x2)=yc_all(:,:,ind);
cfrac_all2(y1:y2,x1:x2,:)=cfrac_all(:,:,:,ind);rA_all2(y1:y2,x1:x2)=rA_all(:,:,ind);
bathy_all2(y1:y2,x1:x2)=bathy_all(:,:,ind);

x1=x2+1; x2=x2+txd;
end

y1=y2+1; y2=y2+tyd;
end
% --------------------------------- % 

% Redifine variables
cfrac=cfrac_all2; clear cfrac_all cfrac_all2; 
rA=rA_all2; clear rA_all; rA_all2;
xc=xc_all2; clear xc_all xc_all2;
yc=yc_all2; clear yc_all yc_all2; 
bathy=bathy_all2; clear bathy_all bathy_all2;
bathy(find(bathy==0))=NaN; % replace land values with NaN

% Define masks
mask3d=cfrac; mask3d(find(cfrac==0))=NaN; mask3d(find(~isnan(mask3d)))=1;
mask2d=mask3d(:,:,1);

% Calculate Cell Volume 
vol=[];
for zz=1:size(cfrac,3)
vol(:,:,zz)=rA*depth_dz(zz).*cfrac(:,:,zz);
end

% =================================================================== % 

% Define Model Field Dimensions
xdim=size(bathy,2);
ydim=size(bathy,1);
zdim=length(depth_uv);
tdim=244;

% Model Start Time 
msd=datenum(1992,12,31);	% conceptually midnight (first instantaneous output would be 01/01/1993 midnight)

display('Pause');pause

% ###################################################################### %
% ######################## READ 3D Monthly Ouput ####################### %
% ###################################################################### %

lp_npp_field_all=[];
out_time=[];

% ******** Time Record Loop ******* % 
for tt=1:244
display(['Record = ',num2str(tt)]);tic

lp_npp_field=[];

% ****** Directory Loop ********** % 
for dd=1:length(dir_names)
files_bio3D=dir([path1,dir_names(dd).name,'/*3D_bio*']);
nc_file=[path1,dir_names(dd).name,'/',files_bio3D(1).name];
ncid=netcdf.open(nc_file,'NOWRITE');
varid=netcdf.inqVarID(ncid,'lp_NPP_f');lp_npp=netcdf.getVar(ncid,varid,[0,0,0,tt-1],[txd,tyd,tzd,1]);lp_npp=double(lp_npp);
varid=netcdf.inqVarID(ncid,'T');rec=netcdf.getVar(ncid,varid,[tt-1],[1]);rec=double(rec)/86400+msd;
netcdf.close(ncid);

% --------- Rotate ------- % 
lp_npp2=[];
for zz=1:zdim
lp_npp2(:,:,zz)=lp_npp(:,:,zz)';
end
lp_npp=lp_npp2; clear npp2;
% ----------------------------- % 

% Save Tile 
lp_npp_field(:,:,:,dd)=lp_npp;
end 	% Directory Loop 
% ***************************************** % 

% ---------- Organize into single array ---------- % 
lp_npp_field2=[];

y1=1; y2=tyd;
for yy=1:num_ty
x1=1; x2=txd;
for xx=1:num_tx
ind=find(xlim_all(:,1)==xlim(xx,1) & xlim_all(:,2)==xlim(xx,2)...
& ylim_all(:,1)==ylim(yy,1) & ylim_all(:,2)==ylim(yy,2));
lp_npp_field2(y1:y2,x1:x2,:)=lp_npp_field(:,:,:,ind);
x1=x2+1; x2=x2+txd;
end
y1=y2+1; y2=y2+tyd;
end
toc

% Redifine Variables and Apply Mask 
lp_npp_field=lp_npp_field2.*mask3d; clear lp_npp_field2;

% Apply scaling (x48) to daily units
lp_npp_field=lp_npp_field*48;

% Save field
lp_npp_field_all(:,:,:,tt)=lp_npp_field;
out_time(1,tt)=rec;
end 

% ####################################################################### % 
% ####################################################################### % 
% #######################################################################
display('Done Reading Fields');

% Get correct time stamp
time=[msd, out_time];
time_final=[];
for ii=1:length(time)-1
time_final(1,ii)=(time(ii)+time(ii+1))/2;
end 

% Get January fields 
january_lp_npp=lp_npp_field_all(:,:,:,1); date_all=[];
count=1;
for yy=1993:2012
target_date=datenum(yy,01,15);
time_diff=time_final-target_date;

ind=find(time_diff==0);
if isempty(ind)
t1=max(time_final(find(time_diff<0)));
t2=min(time_final(find(time_diff>0)));
ind=find(time_final==t1); field1=lp_npp_field_all(:,:,:,ind);
ind=find(time_final==t2); field2=lp_npp_field_all(:,:,:,ind);
lp_npp=field1+(target_date-t1).*(field2-field1)./(t2-t1);
display('Interpolation');
else 
lp_npp=lp_npp_field_all(:,:,:,ind);
display('Match');
end 

january_lp_npp(:,:,:,count)=lp_npp;
date_all(1,count)=target_date;
count=count+1;
end

% Plot
figure
for ii=1:size(january_lp_npp,4)
pcolor(xc,yc,january_lp_npp(:,:,1,ii));shading flat;colorbar; caxis([0 0.1])
title(['Date = ',datestr(date_all(ii))])
pause;
clf
end 








