clear all; close all;

%===========================================================
% OBJECTIVE
%   Regrid lat/lon soluble Fe deposition flux 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. Original solfe 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/solfe/')               % soluble Fe data

% load original data on a lat/lon grid
ncload scenario3.nc

%----------------------------------------------------------
% longitude & latitude
%----------------------------------------------------------

% convert data to double
lon       = double(lon);  % [  0     , 358.1250]     
lat       = double(lat);  % [-88.5419,  88.5419]     

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

%----------------------------------------------------------
% monthly soluble iron deposition
%----------------------------------------------------------

% convert data to double 
dataField_latlon = double(hp_cp_dst_solfe_dep_current);   % kg/m2/s

% convert kg/m2/s to molFe/m2/s
dataField_latlon = (1000./55.845) * dataField_latlon; % dim [12 93 192]

% replace missing values with NaNs
dataField_latlon( dataField_latlon == -9999 ) = NaN;

% rearrange dim to be [lon lat month] 
dataField_latlon = permute( dataField_latlon, [3 2 1]);


%----------------------------------------------------------
% regrid 2D soluble iron dep from lat/lon to LLC90 one
% month at a time
%----------------------------------------------------------

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

% loop over month
for j = 1:12;
   eval([ 'dataField_llc90_', num2str(j),...
          '= gcmfaces_remap_2d( lon, lat, squeeze(dataField_latlon(:,:,j)), 2);'  ]);
end

% write to binary
write2file( 'mahowald_feDep_llc90_01.bin', convert2gcmfaces( dataField_llc90_1 ) );
write2file( 'mahowald_feDep_llc90_02.bin', convert2gcmfaces( dataField_llc90_2 ) );
write2file( 'mahowald_feDep_llc90_03.bin', convert2gcmfaces( dataField_llc90_3 ) );
write2file( 'mahowald_feDep_llc90_04.bin', convert2gcmfaces( dataField_llc90_4 ) );
write2file( 'mahowald_feDep_llc90_05.bin', convert2gcmfaces( dataField_llc90_5 ) );
write2file( 'mahowald_feDep_llc90_06.bin', convert2gcmfaces( dataField_llc90_6 ) );
write2file( 'mahowald_feDep_llc90_07.bin', convert2gcmfaces( dataField_llc90_7 ) );
write2file( 'mahowald_feDep_llc90_08.bin', convert2gcmfaces( dataField_llc90_8 ) );
write2file( 'mahowald_feDep_llc90_09.bin', convert2gcmfaces( dataField_llc90_9 ) );
write2file( 'mahowald_feDep_llc90_10.bin', convert2gcmfaces( dataField_llc90_10 ) );
write2file( 'mahowald_feDep_llc90_11.bin', convert2gcmfaces( dataField_llc90_11 ) );
write2file( 'mahowald_feDep_llc90_12.bin', convert2gcmfaces( dataField_llc90_12 ) );

%----------------------------------------------------------
% visual check: plot January soluble Fe deposition (molFe/m2/s)
%----------------------------------------------------------
face1 = dataField_llc90_1.f1;
face2 = dataField_llc90_1.f2;
face3 = dataField_llc90_1.f3;
face4 = dataField_llc90_1.f4;
face5 = dataField_llc90_1.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;

return


