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 |
|
|
% 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 |
|
|
|
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 |
|
|
FocnNew = mean(dataCplSelect.([mncfld{ifield},'tave']),3); |
66 |
|
|
FocnNew = DataCorrections(FocnNew,mncfld{ifield},gravity,RAC); |
67 |
|
|
if FocnOnly |
68 |
|
|
Focn = FocnNew; |
69 |
|
|
else |
70 |
|
|
Focn = (gWeight*FocnOld + FocnNew)/(gWeight+1); |
71 |
|
|
end |
72 |
|
|
fid=fopen([outDir,'/',FocnRoot,'.',fields{ifield},... |
73 |
|
|
'.bin.',num2str(period)],'w','b'); |
74 |
|
|
fwrite(fid,Focn,'real*8'); |
75 |
|
|
fclose(fid); |
76 |
|
|
dataFocn.(mncfld{ifield}) = Focn; |
77 |
|
|
end |
78 |
|
|
|
79 |
|
|
|
80 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
81 |
|
|
% Prepare focn as needed % |
82 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
83 |
|
|
|
84 |
|
|
if ~FocnOnly |
85 |
|
|
|
86 |
|
|
years = [iitC:ditC:fitC].*dtC./31104000; |
87 |
|
|
yrfrc = rem(years,floor(years)); |
88 |
|
|
yrfrc(yrfrc == 0) = 1; |
89 |
|
|
uyrfrc = sort(unique(yrfrc)); |
90 |
|
|
length(yrfrc) |
91 |
|
|
length(uyrfrc) |
92 |
|
|
|
93 |
|
|
for ifield = 1:length(mncfld) |
94 |
|
|
|
95 |
|
|
% Make yearly forcing file. |
96 |
|
|
focnFull = dataCplSelect.([mncfld{ifield},'tave']); |
97 |
|
|
focn = NaN.*zeros([size(RAC,1),size(RAC,2),length(uyrfrc)]); |
98 |
|
|
for itime = 1:length(uyrfrc) |
99 |
|
|
index = find(uyrfrc(itime) == yrfrc); |
100 |
|
|
focn(:,:,itime) = mean(focnFull(:,:,index),3); |
101 |
|
|
end |
102 |
|
|
if ~isempty(find(isnan(focn(:)))) |
103 |
|
|
error('NaNs in ocean forcing files!!'); |
104 |
|
|
end |
105 |
|
|
|
106 |
|
|
% Remove mean, add Focn |
107 |
|
|
SubMean = mean(focn,3); |
108 |
|
|
AddMean = dataFocn.(mncfld{ifield}); |
109 |
|
|
for itime = 1:size(focn,3) |
110 |
|
|
focn(:,:,itime) = focn(:,:,itime) - SubMean + AddMean; |
111 |
|
|
end |
112 |
|
|
|
113 |
|
|
% Where ice, make Tau_x, Tau_y equal zero. |
114 |
|
|
if isequal(mncfld{ifield},'SIC') |
115 |
|
|
iceMask = (focn==0); |
116 |
|
|
elseif ismember(mncfld{ifield},{'TX','TY'}) |
117 |
|
|
focn = focn.*iceMask; |
118 |
|
|
end |
119 |
|
|
|
120 |
|
|
|
121 |
|
|
focnmean = mean(focn,3); |
122 |
|
|
focnmean = sum(focnmean(:).*RAC(:))./sum(RAC(:)); |
123 |
|
|
disp([mncfld{ifield},' mean: ',num2str(focnmean)]); |
124 |
|
|
|
125 |
|
|
fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'w','b'); |
126 |
|
|
fwrite(fid,focn,'real*8'); |
127 |
|
|
fclose(fid); |
128 |
|
|
end |
129 |
|
|
end |
130 |
|
|
|
131 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
132 |
|
|
% Inspect Focn % |
133 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
134 |
|
|
|
135 |
|
|
if InspectFocn |
136 |
|
|
for ifield = 1:length(mncfld) |
137 |
|
|
fid=fopen([outDir,'/',FocnRoot,'.',fields{ifield},... |
138 |
|
|
'.bin.',num2str(period)],'r','b'); |
139 |
|
|
Focn = reshape(fread(fid,'real*8'),size(RAC)); |
140 |
|
|
fclose(fid); |
141 |
|
|
figure(ifield); |
142 |
|
|
merccube(XG,YG,Focn); colorbar; |
143 |
|
|
title(mncfld{ifield}); |
144 |
|
|
end |
145 |
|
|
end |
146 |
|
|
|
147 |
|
|
if Inspectfocn |
148 |
|
|
for ifield = 1:length(mncfld) |
149 |
|
|
figure(ifield); orient tall; |
150 |
|
|
fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'r','b'); |
151 |
|
|
focn = reshape(fread(fid,'real*8'),[size(RAC),12]); |
152 |
|
|
fclose(fid); |
153 |
|
|
crange = [min(focn(:)),max(focn(:))]; |
154 |
|
|
for ii = 1:size(focn,3) |
155 |
|
|
subplot(4,3,ii); merccube(XG,YG,focn(:,:,ii)); |
156 |
|
|
shading flat; caxis(crange); colorbar; |
157 |
|
|
title([mncfld{ifield},' month ',num2str(ii)]); |
158 |
|
|
end |
159 |
|
|
end |
160 |
|
|
end |
161 |
|
|
|
162 |
|
|
|
163 |
|
|
% Field dependent data corrections. |
164 |
|
|
function CorrectedData = ... |
165 |
|
|
DataCorrections(DataToCorrect,CorrectionField,gravity,RAC) |
166 |
|
|
if isequal(CorrectionField,'SIC') |
167 |
|
|
CorrectedData = gravity.*DataToCorrect; |
168 |
|
|
elseif isequal(CorrectionField,'FW') |
169 |
|
|
offset = sum(DataToCorrect(:).*RAC(:))./sum(RAC(:)); |
170 |
|
|
CorrectedData = DataToCorrect - offset; |
171 |
|
|
else |
172 |
|
|
CorrectedData = DataToCorrect; |
173 |
|
|
end |
174 |
|
|
|