enderton/PeriodicCoupling/CalcOcnForcingFields.m

function CalcOcnForcingFields(cplFiles,gridFiles,prevFocnFiles,period,...
                              outDir,FocnRoot,focnRoot,iitC,fitC,ditC,...
                              dtC,gWeight,FocnOnly,InspectFocn,Inspectfocn)

fields = {'icL','taux','tauy','EmP','Qnet'};
mncfld = {'SIC','TX'  ,'TY'  ,'FW' ,'HF'  };

gravity = 9.81;

warning off MATLAB:MKDIR:DirectoryExists;
mkdir(outDir);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 Load Data                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

datacpl = rdmnc(cplFiles,'iter','HFtave','TXtave',...
                'TYtave','FWtave','SICtave','T');
datagrd = rdmnc(gridFiles,'rA','XG','YG','XC','YC','HFacC');
RAC = datagrd.rA;
HFacC = datagrd.HFacC;
XC = datagrd.XC; XG = datagrd.XG(1:size(RAC,1),1:size(RAC,2));
YC = datagrd.YC; YG = datagrd.YG(1:size(RAC,1),1:size(RAC,2));
% save('DataPerCpl.mat','datacpl','RAC','XG','YG','XC','YC','HFacC');
% load('DataPerCpl.mat');
HFacC_Top = HFacC(:,:,1);

% If a new focn is to be calculated, the old Focn must be loaded.
if ~FocnOnly
    for ifield=1:length(mncfld)
        filename = [prevFocnFiles,'/',FocnRoot,'.',...
                   fields{ifield},'.bin.',num2str(period-1)];
        fid=fopen(filename,'r','b');
        dataFocnOld.(mncfld{ifield}) = ...
            reshape(fread(fid,'real*8'),size(RAC));
        fclose(fid);
    end
end
% merccube_mod(XG,YG,dataFocnOld.HF);
% aaa = dataFocnOld.HF;
% sum(aaa(:).*RAC(:))./sum(RAC(:))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      Select Coupled Field Data                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

test  = ismember([iitC:ditC:fitC],datacpl.iter);
index = ismember(datacpl.iter,[iitC:ditC:fitC]);
if ~isempty(find(test==0))
    error('Specified indecies not present in cpl_tave.* files!!');
end
for ifield = 1:length(mncfld)
    dataCplSelect.([mncfld{ifield},'tave']) = ...
        datacpl.([mncfld{ifield},'tave'])(:,:,index);
end
% merccube_mod(XG,YG,dataCplSelect.HFtave(:,:,end));
% aaa = dataCplSelect.HFtave(:,:,end);
% sum(aaa(:).*RAC(:))./sum(RAC(:))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Prepare Focn                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% What if it is not in year increments?  Need a more sophisticated
% treatment of finding Focn.  Should compute monthly mean before.  Test
% with various sorts of coupled field outputs.
for ifield = 1:length(mncfld)
    FocnNew = HFacC_Top.*mean(dataCplSelect.([mncfld{ifield},'tave']),3);
    FocnNew = DataCorrections(FocnNew,mncfld{ifield},gravity,RAC,HFacC_Top);
    if FocnOnly
        Focn = FocnNew;
    else
        FocnOld = HFacC_Top.*dataFocnOld.(mncfld{ifield});
        if isequal(mncfld{ifield},'FW')
            FocnOld = DataCorrections(FocnOld,mncfld{ifield},...
                                      gravity,RAC,HFacC_Top);
        end
        Focn = (gWeight*FocnOld + FocnNew)/(gWeight+1);
    end
    fid=fopen([outDir,'/',FocnRoot,'.',fields{ifield},...
               '.bin.',num2str(period)],'w','b');
    fwrite(fid,Focn,'real*8');
    fclose(fid);
    dataFocn.(mncfld{ifield}) = Focn;
end
% merccube_mod(XG,YG,dataFocn.HF);
% aaa = dataFocn.HF(:,:,end);
% sum(aaa(:).*RAC(:))./sum(RAC(:))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         Prepare focn as needed                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Must account for land!!!!
if ~FocnOnly
    
    years = [iitC:ditC:fitC].*dtC./31104000;
    yrfrc = rem(years,floor(years));
    yrfrc(yrfrc == 0) = 1;
    uyrfrc = sort(unique(yrfrc));
    disp(['Cpl. data instances:  ',num2str(length(yrfrc))]);
    disp(['Cpl. data per year forcing:  ',length(uyrfrc)]);
    
    for ifield = 1:length(mncfld)
        
        % Make yearly forcing file (usually twelve months).
        focnFull = dataCplSelect.([mncfld{ifield},'tave']);
        focn = NaN.*zeros([size(RAC,1),size(RAC,2),length(uyrfrc)]);
        for itime = 1:length(uyrfrc)
            index = find(uyrfrc(itime) == yrfrc);
            focn(:,:,itime) = mean(focnFull(:,:,index),3);
        end
        if ~isempty(find(isnan(focn(:))))
            error('NaNs in ocean forcing files!!');
        end
        
        % Account for bathymetry.
        for itime = 1:size(focn,3)
            focn(:,:,itime) = HFacC_Top.*focn(:,:,itime);
        end
        
        % Remove time mean from focn, add Focn
        SubMean = mean(focn,3);
        AddMean = dataFocn.(mncfld{ifield});
        for itime = 1:size(focn,3)
            focn(:,:,itime) = focn(:,:,itime) - SubMean + AddMean;
        end
        
        % Where ice, make Tau_x, Tau_y equal zero.
        if isequal(mncfld{ifield},'SIC')
            iceMask = (focn==0);
        elseif ismember(mncfld{ifield},{'TX','TY'})
            focn = focn.*iceMask;
        end
        
        focnmean = mean(focn,3);
        focnmean = sum(focnmean(:).*RAC(:))./sum(RAC(:));
        disp([mncfld{ifield},' mean:  ',num2str(focnmean)]);
        
        fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'w','b');
        fwrite(fid,focn,'real*8');
        fclose(fid);
        
        aaa = mean(focn,3);
        disp(['Number of land cells that do not have zero values:  ',...
              num2str(length(find(aaa(find(~HFacC_Top))~=0)))]);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Inspect Focn                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if InspectFocn
    for ifield = 1:length(mncfld)
        fid=fopen([outDir,'/',FocnRoot,'.',fields{ifield},...
                   '.bin.',num2str(period)],'r','b');
        Focn = reshape(fread(fid,'real*8'),size(RAC));
        fclose(fid);
        figure;
        merccube_mod(XG,YG,Focn); colorbar;
        title(mncfld{ifield});
    end
end

if Inspectfocn
    for ifield = 1:length(mncfld)
        figure; orient tall;
        fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'r','b');
        focn = reshape(fread(fid,'real*8'),[size(RAC),12]);
        fclose(fid);
        crange = [min(focn(:)),max(focn(:))];
        for ii = 1:size(focn,3)
            subplot(4,3,ii); merccube_mod(XG,YG,focn(:,:,ii));
            shading flat; caxis(crange); colorbar;
            title([mncfld{ifield},' month ',num2str(ii)]);
        end
    end
end


% Field dependent data corrections.
function CorrectedData = ...
  DataCorrections(DataToCorrect,CorrectionField,gravity,RAC,HFacC_Top)
    index = find(HFacC_Top);
    correctdim = size(DataToCorrect);
    DataToCorrect = DataToCorrect(:);
    if     isequal(CorrectionField,'SIC')
        CorrectedData(index) = gravity.*DataToCorrect(index);
    elseif isequal(CorrectionField,'FW')
        DataToCorrect(index);
        offset = sum(DataToCorrect(index).*RAC(index))./sum(RAC(index));
        CorrectedData(index) = DataToCorrect(index) - offset;
    else
        CorrectedData(index) = DataToCorrect(index);
    end
    CorrectedData = reshape(CorrectedData,correctdim);
        
1	enderton	1.1	function CalcOcnForcingFields(cplFiles,gridFiles,prevFocnFiles,period,...
2			outDir,FocnRoot,focnRoot,iitC,fitC,ditC,...
3			dtC,gWeight,FocnOnly,InspectFocn,Inspectfocn)
4
5			fields = {'icL','taux','tauy','EmP','Qnet'};
6			mncfld = {'SIC','TX' ,'TY' ,'FW' ,'HF' };
7
8			gravity = 9.81;
9
10			warning off MATLAB:MKDIR:DirectoryExists;
11			mkdir(outDir);
12
13			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14			% Load Data %
15			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16
17	enderton	1.2	datacpl = rdmnc(cplFiles,'iter','HFtave','TXtave',...
18	enderton	1.4	'TYtave','FWtave','SICtave','T');
19			datagrd = rdmnc(gridFiles,'rA','XG','YG','XC','YC','HFacC');
20	enderton	1.2	RAC = datagrd.rA;
21	enderton	1.4	HFacC = datagrd.HFacC;
22	enderton	1.2	XC = datagrd.XC; XG = datagrd.XG(1:size(RAC,1),1:size(RAC,2));
23			YC = datagrd.YC; YG = datagrd.YG(1:size(RAC,1),1:size(RAC,2));
24	enderton	1.4	% save('DataPerCpl.mat','datacpl','RAC','XG','YG','XC','YC','HFacC');
25			% load('DataPerCpl.mat');
26			HFacC_Top = HFacC(:,:,1);
27	enderton	1.1
28			% If a new focn is to be calculated, the old Focn must be loaded.
29			if ~FocnOnly
30			for ifield=1:length(mncfld)
31			filename = [prevFocnFiles,'/',FocnRoot,'.',...
32			fields{ifield},'.bin.',num2str(period-1)];
33			fid=fopen(filename,'r','b');
34			dataFocnOld.(mncfld{ifield}) = ...
35			reshape(fread(fid,'real*8'),size(RAC));
36			fclose(fid);
37			end
38			end
39	enderton	1.4	% merccube_mod(XG,YG,dataFocnOld.HF);
40			% aaa = dataFocnOld.HF;
41			% sum(aaa(:).*RAC(:))./sum(RAC(:))
42	enderton	1.1
43
44			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45			% Select Coupled Field Data %
46			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47
48			test = ismember([iitC:ditC:fitC],datacpl.iter);
49			index = ismember(datacpl.iter,[iitC:ditC:fitC]);
50			if ~isempty(find(test==0))
51			error('Specified indecies not present in cpl_tave.* files!!');
52			end
53			for ifield = 1:length(mncfld)
54			dataCplSelect.([mncfld{ifield},'tave']) = ...
55			datacpl.([mncfld{ifield},'tave'])(:,:,index);
56			end
57	enderton	1.4	% merccube_mod(XG,YG,dataCplSelect.HFtave(:,:,end));
58			% aaa = dataCplSelect.HFtave(:,:,end);
59			% sum(aaa(:).*RAC(:))./sum(RAC(:))
60	enderton	1.1
61
62			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63			% Prepare Focn %
64			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65
66			% What if it is not in year increments? Need a more sophisticated
67			% treatment of finding Focn. Should compute monthly mean before. Test
68			% with various sorts of coupled field outputs.
69			for ifield = 1:length(mncfld)
70	enderton	1.4	FocnNew = HFacC_Top.*mean(dataCplSelect.([mncfld{ifield},'tave']),3);
71			FocnNew = DataCorrections(FocnNew,mncfld{ifield},gravity,RAC,HFacC_Top);
72	enderton	1.1	if FocnOnly
73			Focn = FocnNew;
74			else
75	enderton	1.4	FocnOld = HFacC_Top.*dataFocnOld.(mncfld{ifield});
76	enderton	1.3	if isequal(mncfld{ifield},'FW')
77	enderton	1.4	FocnOld = DataCorrections(FocnOld,mncfld{ifield},...
78			gravity,RAC,HFacC_Top);
79	enderton	1.3	end
80	enderton	1.1	Focn = (gWeight*FocnOld + FocnNew)/(gWeight+1);
81			end
82			fid=fopen([outDir,'/',FocnRoot,'.',fields{ifield},...
83			'.bin.',num2str(period)],'w','b');
84			fwrite(fid,Focn,'real*8');
85			fclose(fid);
86			dataFocn.(mncfld{ifield}) = Focn;
87			end
88	enderton	1.4	% merccube_mod(XG,YG,dataFocn.HF);
89			% aaa = dataFocn.HF(:,:,end);
90			% sum(aaa(:).*RAC(:))./sum(RAC(:))
91	enderton	1.1
92
93			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94			% Prepare focn as needed %
95			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96
97	enderton	1.4	% Must account for land!!!!
98	enderton	1.1	if ~FocnOnly
99
100			years = [iitC:ditC:fitC].*dtC./31104000;
101			yrfrc = rem(years,floor(years));
102			yrfrc(yrfrc == 0) = 1;
103			uyrfrc = sort(unique(yrfrc));
104	enderton	1.4	disp(['Cpl. data instances: ',num2str(length(yrfrc))]);
105			disp(['Cpl. data per year forcing: ',length(uyrfrc)]);
106	enderton	1.1
107			for ifield = 1:length(mncfld)
108
109	enderton	1.4	% Make yearly forcing file (usually twelve months).
110	enderton	1.1	focnFull = dataCplSelect.([mncfld{ifield},'tave']);
111			focn = NaN.*zeros([size(RAC,1),size(RAC,2),length(uyrfrc)]);
112			for itime = 1:length(uyrfrc)
113			index = find(uyrfrc(itime) == yrfrc);
114			focn(:,:,itime) = mean(focnFull(:,:,index),3);
115			end
116			if ~isempty(find(isnan(focn(:))))
117			error('NaNs in ocean forcing files!!');
118			end
119
120	enderton	1.4	% Account for bathymetry.
121			for itime = 1:size(focn,3)
122			focn(:,:,itime) = HFacC_Top.*focn(:,:,itime);
123			end
124
125			% Remove time mean from focn, add Focn
126	enderton	1.1	SubMean = mean(focn,3);
127			AddMean = dataFocn.(mncfld{ifield});
128			for itime = 1:size(focn,3)
129			focn(:,:,itime) = focn(:,:,itime) - SubMean + AddMean;
130			end
131
132			% Where ice, make Tau_x, Tau_y equal zero.
133			if isequal(mncfld{ifield},'SIC')
134			iceMask = (focn==0);
135			elseif ismember(mncfld{ifield},{'TX','TY'})
136			focn = focn.*iceMask;
137			end
138
139			focnmean = mean(focn,3);
140			focnmean = sum(focnmean(:).*RAC(:))./sum(RAC(:));
141			disp([mncfld{ifield},' mean: ',num2str(focnmean)]);
142
143			fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'w','b');
144			fwrite(fid,focn,'real*8');
145			fclose(fid);
146	enderton	1.4
147			aaa = mean(focn,3);
148			disp(['Number of land cells that do not have zero values: ',...
149			num2str(length(find(aaa(find(~HFacC_Top))~=0)))]);
150	enderton	1.1	end
151			end
152
153	enderton	1.4
154	enderton	1.1	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155			% Inspect Focn %
156			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
157
158			if InspectFocn
159			for ifield = 1:length(mncfld)
160			fid=fopen([outDir,'/',FocnRoot,'.',fields{ifield},...
161			'.bin.',num2str(period)],'r','b');
162			Focn = reshape(fread(fid,'real*8'),size(RAC));
163			fclose(fid);
164	enderton	1.4	figure;
165			merccube_mod(XG,YG,Focn); colorbar;
166	enderton	1.1	title(mncfld{ifield});
167			end
168			end
169
170			if Inspectfocn
171			for ifield = 1:length(mncfld)
172	enderton	1.4	figure; orient tall;
173	enderton	1.1	fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'r','b');
174			focn = reshape(fread(fid,'real*8'),[size(RAC),12]);
175			fclose(fid);
176			crange = [min(focn(:)),max(focn(:))];
177			for ii = 1:size(focn,3)
178	enderton	1.4	subplot(4,3,ii); merccube_mod(XG,YG,focn(:,:,ii));
179	enderton	1.1	shading flat; caxis(crange); colorbar;
180			title([mncfld{ifield},' month ',num2str(ii)]);
181			end
182			end
183			end
184
185
186			% Field dependent data corrections.
187			function CorrectedData = ...
188	enderton	1.4	DataCorrections(DataToCorrect,CorrectionField,gravity,RAC,HFacC_Top)
189			index = find(HFacC_Top);
190			correctdim = size(DataToCorrect);
191			DataToCorrect = DataToCorrect(:);
192	enderton	1.1	if isequal(CorrectionField,'SIC')
193	enderton	1.4	CorrectedData(index) = gravity.*DataToCorrect(index);
194	enderton	1.1	elseif isequal(CorrectionField,'FW')
195	enderton	1.4	DataToCorrect(index);
196			offset = sum(DataToCorrect(index).*RAC(index))./sum(RAC(index));
197			CorrectedData(index) = DataToCorrect(index) - offset;
198	enderton	1.1	else
199	enderton	1.4	CorrectedData(index) = DataToCorrect(index);
200	enderton	1.1	end
201	enderton	1.4	CorrectedData = reshape(CorrectedData,correctdim);
202	enderton	1.2