clear all; close all;

%===========================================================
% OBJECTIVE
%   Regrid lat/lon 3D total iron from IGSM model runs 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
%   01 Apr 2015 - H. Amos - v 1.0 created. IGSM model data from
%                           
% 
% 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/') )       % for gcmfaces package
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( genpath('~/MITgcm/darwin_data/from_StephD/') ) % soluble IGSM data

% load IGSM grid information
ncload grid_igsm.nc

% convert grid information to double
lat   = double(Y); % dim [1 90 ]
lon   = double(X); % dim [1 144], range [1.25 358.75]
depth = double(Z); % dim [1 22 ]

% put lat/lon on a grid (i.e., 2D matrix)
[lat, lon] = meshgrid(lat, lon);

% load 3D IGSM data field
dataField_latlon = read2memory('~/MITgcm/darwin_data/from_StephD/IGSM/fet_y1990.data',...
                               [144 90 22]);

% make sure data is double
dataField_latlon = double( dataField_latlon );

% replacing missing values (0s) with NaNs
dataField_latlon( dataField_latlon == 0 ) = NaN;

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

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

% write LLC90 gridded data to binary for input to MITgcm
write2file( 'igsm_fet_llc90.bin' , convert2gcmfaces( dataField_llc90 ) );


disp('Success!')
return
