clear all; close all;

%===========================================================
% OBJECTIVE
%   Regrid lat/lon SeaWIFS data to the ECCO v4 LLC90 grid. 
%
% NOTES
%   (01) This routine requires the following external MATLAB
%        files:
%          - rdmds.m
%          - gcmfaces/ package from Gael Forget (MIT)
%
% REVISION HISTORY
%   30 Mar 2015 - H. Amos - v 1.0 created. Original SeaWIFS code
%                           written by Yanxu Zhang (Harvard)
%                           and then I've adapted it based on 
%                           woa13_regrid.m from Gael Forget (MIT).
% 
% Helen M. Amos
% hamos@hsph.harvard.edu
% Harvard T. H. Chan School of Public Health
%===========================================================

% add paths 
addpath('~/MATLAB')                                  % for rdms.m
addpath('~/MATLAB/gcmfaces/gcmfaces/')               % for gcmfaces package
addpath( genpath('~/MATLAB/gcmfaces/gcmfaces/') )
addpath('~/MATLAB/gcmfaces/gcmfaces/gcmfaces_calc/') % for gcmfaces_remap_2d
addpath('~/MATLAB/gcmfaces/gcmfaces/gcmfaces_IO/')   % for grid_load.m
addpath('~/MATLAB/gcmfaces/MITprof/profiles_misc/')  % for ncload
addpath('~/MITgcm/darwin_data/plankton/tempdata/')   % SeaWIFS Chl data

% reformat SeaWIFS lat/lon for LCC90 regridding function
lon        = linspace(1,360,360)-1-179.5; % make a row vector of longitudes
lon        = lon'*ones(1,180);       % vector to matrix
lon(lon<0) = lon(lon<0)+360;         % convert [-180 180] to [0 360]
lon        = circshift(lon,[180 0]); % shift position of vector elements by 180

lat =linspace(1,180,180)-1-89.5;     % make a row vector of latitudes
lat        = ones(360,1)*lat;        % vector to marix

%-----------------------------------------------------------
% load 2D SeaWIFS Chl data, then reformat
%-----------------------------------------------------------

% load ascii file pre-processed by Yanxu Zhang
load chl.txt;                             

% rotate matrix
chl_latlon_2d                      = chl';    

% replace missing numbers with NaN                
chl_latlon_2d(chl_latlon_2d == 0.) = NaN; 

% transform 2D to 3D by copying surface layer
chl_latlon_3d                      = repmat( chl_latlon_2d, [1 1 50]); 

% shift position of array elements by 180
chl_latlon_3d                      = circshift(chl_latlon_3d, [180 0 0]);  

% create a dummy depth row vector of depth (m)
depth = linspace(0,4500,50);

%-----------------------------------------------------------
% regrid lat/lon to LLC90
%-----------------------------------------------------------

% LLC90 grid information
grid_load('~/MATLAB/gcmfaces/GRID/',5,'compact');

% regrid 3D lat/lon SeaWIFS data to LLC90 
dataField_llc90 = gcmfaces_remap_3d( lon, lat, depth, chl_latlon_3d );
 

% save to .mat, for debugging (hma, 19 mar 2015)
save('TEST_seawifs.mat','dataField_llc90')

% 1. Convert woa_poa_llc90 to a gcmfaces object.
% 2. Write LLC90 gridded WOA data to binary for input to MITgcm. 
write2file( 'seawifs_plankton_llc90.bin' , convert2gcmfaces(dataField_llc90) );


%-----------------------------------------------------------
% visual check: plot surface layer of 3D field
%-----------------------------------------------------------
face1 = dataField_llc90.f1;
face2 = dataField_llc90.f2;
face3 = dataField_llc90.f3;
face4 = dataField_llc90.f4;
face5 = dataField_llc90.f5;

figure(1)
subplot(2,3,1)
pcolor(face1(:,:,1)); shading flat;
subplot(2,3,2)
pcolor(face2(:,:,1)); shading flat;
subplot(2,3,3)
pcolor(face3(:,:,1)); shading flat;
subplot(2,3,4)
pcolor(face4(:,:,1)); shading flat;
subplot(2,3,5)
pcolor(face5(:,:,1)); shading flat;

%-----------------------------------------------------------
% transfer data to have a specific mean value, replace NaNs
% with zeros, then write binary
%-----------------------------------------------------------

faceAbbrev = {'f1','f2','f3','f4','f4'};

% all plankton data have units of mmol/m3

%------------------------------
% phy
%------------------------------
phy_llc90 = dataField_llc90./nanmean( nanmean( dataField_llc90 )  ) * 0.01;

% zero out subsurface layers and replace NaNs with zeros 
% in the surface
for n = 1:numel(faceAbbrev); % loop over faces
  clear temp; % for safety's sake
  eval([ 'temp = phy_llc90.' faceAbbrev{n} ';' ]);
  temp(:,:,2:end)   = 0.;
  temp(isnan(temp)) = 0.;
  eval([ 'phy_llc90.' faceAbbrev{n} '= temp;' ]);
end

% write  binary
write2file( 'phy_seawifs_llc90.bin', convert2gcmfaces( phy_llc90 ));

%------------------------------
% zooc
%------------------------------
zooc_llc90 = dataField_llc90./nanmean( nanmean( dataField_llc90 )  ) * 0.21;

% replace NaNs with zeros in subsurface layers
for n = 1:numel(faceAbbrev); % loop over faces
  clear temp; % for safety's sake
  eval([ 'temp = zooc_llc90.' faceAbbrev{n} ';' ]);
  temp(:,:,2:end)   = 0.;
  temp(isnan(temp)) = 0.;
  eval([ 'zooc_llc90.' faceAbbrev{n} '= temp;' ]);
end

% write  binary
write2file( 'zooc_seawifs_llc90.bin', convert2gcmfaces( zooc_llc90 ));

%------------------------------
% zoofe 
%------------------------------
zoofe_llc90 = dataField_llc90./nanmean( nanmean( dataField_llc90 )  ) * 1.83e-6;

% replace NaNs with zeros in subsurface layers
for n = 1:numel(faceAbbrev); % loop over faces
  clear temp; % for safety's sake
  eval([ 'temp = zoofe_llc90.' faceAbbrev{n} ';' ]);
  temp(:,:,2:end)   = 0.;
  temp(isnan(temp)) = 0.;
  eval([ 'zoofe_llc90.' faceAbbrev{n} '= temp;' ]);
end

% write  binary
write2file( 'zoofe_seawifs_llc90.bin', convert2gcmfaces( zoofe_llc90 ));

%------------------------------
% zoon 0.03
%------------------------------
zoon_llc90 = dataField_llc90./nanmean( nanmean( dataField_llc90 )  ) * 0.03;

% replace NaNs with zeros in subsurface layers
for n = 1:numel(faceAbbrev); % loop over faces
  clear temp; % for safety's sake
  eval([ 'temp = zoon_llc90.' faceAbbrev{n} ';' ]);
  temp(:,:,2:end)   = 0.;
  temp(isnan(temp)) = 0.;
  eval([ 'zoon_llc90.' faceAbbrev{n} '= temp;' ]);
end

% write  binary
write2file( 'zoon_seawifs_llc90.bin', convert2gcmfaces( zoon_llc90 ));

%------------------------------
% zoop 0.002
%------------------------------
zoop_llc90 = dataField_llc90./nanmean( nanmean( dataField_llc90 )  ) * 0.002;

% replace NaNs with zeros in subsurface layers
for n = 1:numel(faceAbbrev); % loop over faces
  clear temp; % for safety's sake
  eval([ 'temp = zoop_llc90.' faceAbbrev{n} ';' ]);
  temp(:,:,2:end)   = 0.;
  temp(isnan(temp)) = 0.;
  eval([ 'zoop_llc90.' faceAbbrev{n} '= temp;' ]);
end

% write  binary
write2file( 'zoop_seawifs_llc90.bin', convert2gcmfaces( zoop_llc90 ));

%------------------------------
% zoosi 0.013
%------------------------------
zoosi_llc90 = dataField_llc90./nanmean( nanmean( dataField_llc90 )  ) * 0.013;

% replace NaNs with zeros in subsurface layers
for n = 1:numel(faceAbbrev); % loop over faces
  clear temp; % for safety's sake
  eval([ 'temp = zoosi_llc90.' faceAbbrev{n} ';' ]);
  temp(:,:,2:end)   = 0.;
  temp(isnan(temp)) = 0.;
  eval([ 'zoosi_llc90.' faceAbbrev{n} '= temp;' ]);
end

% write  binary
write2file( 'zoosi_seawifs_llc90.bin', convert2gcmfaces( zoosi_llc90 ));


disp('Success!')
