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');
datagrd = rdmnc(gridFiles,'rA','XG','YG','XC','YC');
RAC = datagrd.rA;
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');
%load('DataPerCpl.mat');


% 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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            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)
    FocnOld = dataFocnOld.(mncfld{ifield});
    if isequal(mncfld{ifield},'FW')
        FocnOld = DataCorrections(FocnOld,mncfld{ifield},gravity,RAC);
    end
    FocnNew = mean(dataCplSelect.([mncfld{ifield},'tave']),3);
    FocnNew = DataCorrections(FocnNew,mncfld{ifield},gravity,RAC);
    if FocnOnly
        Focn = FocnNew;
    else
        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


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

if ~FocnOnly
    
    years = [iitC:ditC:fitC].*dtC./31104000;
    yrfrc = rem(years,floor(years));
    yrfrc(yrfrc == 0) = 1;
    uyrfrc = sort(unique(yrfrc));
    length(yrfrc)
    length(uyrfrc)
    
    for ifield = 1:length(mncfld)
        
        % Make yearly forcing file.
        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
        
        % Remove mean, 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);
    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(ifield);
        merccube(XG,YG,Focn); colorbar;
        title(mncfld{ifield});
    end
end

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