/[MITgcm]/MITgcm_contrib/enderton/PeriodicCoupling/CalcOcnForcingFields.m
ViewVC logotype

Annotation of /MITgcm_contrib/enderton/PeriodicCoupling/CalcOcnForcingFields.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Wed Oct 5 02:27:21 2005 UTC (19 years, 9 months ago) by enderton
Branch: MAIN
Script for calculating ocean only forcing for ocean only componant of periodic coupling.

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    

  ViewVC Help
Powered by ViewVC 1.1.22