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

Contents of /MITgcm_contrib/enderton/PeriodicCoupling/MakePcForcing.m

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


Revision 1.1 - (show annotations) (download)
Tue Feb 28 23:16:02 2006 UTC (19 years, 4 months ago) by enderton
Branch: MAIN
Modifications to periodic coupling script

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

  ViewVC Help
Powered by ViewVC 1.1.22