%******************************************************************
% This is a pure calculation (no plotting) script: combines  
% modelled thermocline depth anomalies for the first to fifth
% baroclinic vertical modes and convert them to sea level
% anomalies (SLA) during study period of 1961-2002 using
% the coefficients for all vertical mode from Zhu et al (2017)
%
% Input:  mode_[12345]_000000000001_eta_out.nc.trange -from each experiment
%         fit_mm_5.mat -from Zhu et al (2017) 
% Output: eta3d_mm_fit.mat
%	  eta3d_mm_fit.nc
%
% Written by Xiaoting Zhu, 2017 (xzhu@geomar.de)
%******************************************************************

%%%--- BEGIN: global initialization -------------------------
close all;
clear all;
clc;

set(0,'defaultfigurecolor',[1 1 1])
Nmodes = 5;

crfolder = pwd  % current folder: script folder
rtfolder = [crfolder '/..'];  % root folder
dtfolder_supp = [rtfolder '/Greatbatch2018_data/Zhu2017'];
subfolders = {'/NCEP_exp' '/ERA40_exp' '/NCEPzonal_exp' '/ERA40zonal_exp'};

%%%--- Loop: calculate and save for each experiment -------------
for nf = 1:length(subfolders)
 dtfolder = [rtfolder '/Greatbatch2018_data' subfolders{nf}];  % mode...nc.trange folder
 cd(dtfolder)

 %%--- Read model output -----------------------------------------
 lon = ncread('mode_1_000000000001_eta_out.nc.trange', 'LONGITUDE');
 lat = ncread('mode_1_000000000001_eta_out.nc.trange', 'LATITUDE');
 for i = 1:Nmodes
  filename = ['mode_' num2str(i) '_000000000001_eta_out.nc.trange'];
  eval(['mode' num2str(i),'=','ncread(','''',filename,'''',',',...
     '''','eta','''',')',';']); 
 end

 aaa = [];
 ocean = ~isnan(mode2); % ocean mask
 for i = 1:Nmodes
  eval([' aaa =','[aaa,','mode',num2str(i),'(ocean)];']);
  fit_full = aaa;   % only ocean values
 end

 time = ncread('mode_1_000000000001_eta_out.nc.trange', 'TIME');

 %%--- Read coefficients from Zhu et al (2017) ------------------
 fitname = [dtfolder_supp '/fit_mm_' num2str(Nmodes) '.mat'];
 load(fitname);

 %%--- Calculate SLA --------------------------------------------
 eta = NaN(size(mode1));
 mm_slp = nan(size(mode1));
 mm_slp(ocean) = fit_full * coeff;
 eta = mm_slp;
 eta=eta*100; % m->cm

 %%--- Save SLA: both as MAT file and netCDF file ----------------
 save('eta3d_mm_fit.mat','eta','time')
 ncfilename = 'eta3d_mm_fit.nc';
 delete(ncfilename)
 nccreate(ncfilename,'eta','Dimensions',...
                    {'LONGITUDE',length(lon),...
                     'LATITUDE',length(lat),...
                     'TIME',length(time)})
 nccreate(ncfilename,'LONGITUDE','Dimensions',{'LONGITUDE',length(lon)})
 nccreate(ncfilename,'LATITUDE','Dimensions',{'LATITUDE',length(lat)})
 nccreate(ncfilename,'TIME','Dimensions',{'TIME',length(time)})
 ncwriteatt(ncfilename,'LONGITUDE','axis','X')
 ncwriteatt(ncfilename,'LONGITUDE','units','degrees_east')
 ncwriteatt(ncfilename,'LATITUDE','axis','Y')
 ncwriteatt(ncfilename,'LATITUDE','units','degrees_north')
 ncwriteatt(ncfilename,'TIME','calendar','standard')
 ncwriteatt(ncfilename,'TIME','units','seconds since 1957-09-01 00:00:00')
 ncwrite(ncfilename,'eta',eta)
 ncwrite(ncfilename,'LONGITUDE',lon)
 ncwrite(ncfilename,'LATITUDE',lat)
 ncwrite(ncfilename,'TIME',time)
 ncdisp(ncfilename)

end % for nf
