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

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

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

revision 1.2 by enderton, Tue Oct 18 15:51:19 2005 UTC revision 1.4 by enderton, Fri Nov 25 18:28:01 2005 UTC
# Line 15  mkdir(outDir); Line 15  mkdir(outDir);
15  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16    
17  datacpl = rdmnc(cplFiles,'iter','HFtave','TXtave',...  datacpl = rdmnc(cplFiles,'iter','HFtave','TXtave',...
18                           'TYtave','FWtave','SICtave');                  'TYtave','FWtave','SICtave','T');
19  datagrd = rdmnc(gridFiles,'rA','XG','YG','XC','YC');  datagrd = rdmnc(gridFiles,'rA','XG','YG','XC','YC','HFacC');
20  RAC = datagrd.rA;  RAC = datagrd.rA;
21    HFacC = datagrd.HFacC;
22  XC = datagrd.XC; XG = datagrd.XG(1:size(RAC,1),1:size(RAC,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));  YC = datagrd.YC; YG = datagrd.YG(1:size(RAC,1),1:size(RAC,2));
24  %save('DataPerCpl.mat','datacpl','RAC','XG','YG','XC','YC');  % save('DataPerCpl.mat','datacpl','RAC','XG','YG','XC','YC','HFacC');
25  %load('DataPerCpl.mat');  % load('DataPerCpl.mat');
26    HFacC_Top = HFacC(:,:,1);
27    
28  % If a new focn is to be calculated, the old Focn must be loaded.  % If a new focn is to be calculated, the old Focn must be loaded.
29  if ~FocnOnly  if ~FocnOnly
# Line 35  if ~FocnOnly Line 36  if ~FocnOnly
36          fclose(fid);          fclose(fid);
37      end      end
38  end  end
39    % merccube_mod(XG,YG,dataFocnOld.HF);
40    % aaa = dataFocnOld.HF;
41    % sum(aaa(:).*RAC(:))./sum(RAC(:))
42    
43    
44  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 50  for ifield = 1:length(mncfld) Line 54  for ifield = 1:length(mncfld)
54      dataCplSelect.([mncfld{ifield},'tave']) = ...      dataCplSelect.([mncfld{ifield},'tave']) = ...
55          datacpl.([mncfld{ifield},'tave'])(:,:,index);          datacpl.([mncfld{ifield},'tave'])(:,:,index);
56  end  end
57    % merccube_mod(XG,YG,dataCplSelect.HFtave(:,:,end));
58    % aaa = dataCplSelect.HFtave(:,:,end);
59    % sum(aaa(:).*RAC(:))./sum(RAC(:))
60    
61    
62  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 59  end Line 66  end
66  % What if it is not in year increments?  Need a more sophisticated  % What if it is not in year increments?  Need a more sophisticated
67  % treatment of finding Focn.  Should compute monthly mean before.  Test  % treatment of finding Focn.  Should compute monthly mean before.  Test
68  % with various sorts of coupled field outputs.  % with various sorts of coupled field outputs.
   
69  for ifield = 1:length(mncfld)  for ifield = 1:length(mncfld)
70      FocnOld = dataFocnOld.(mncfld{ifield});      FocnNew = HFacC_Top.*mean(dataCplSelect.([mncfld{ifield},'tave']),3);
71      if isequal(mncfld{ifield},'FW')      FocnNew = DataCorrections(FocnNew,mncfld{ifield},gravity,RAC,HFacC_Top);
         FocnOld = DataCorrections(FocnOld,mncfld{ifield},gravity,RAC);  
     end  
     FocnNew = mean(dataCplSelect.([mncfld{ifield},'tave']),3);  
     FocnNew = DataCorrections(FocnNew,mncfld{ifield},gravity,RAC);  
72      if FocnOnly      if FocnOnly
73          Focn = FocnNew;          Focn = FocnNew;
74      else      else
75            FocnOld = HFacC_Top.*dataFocnOld.(mncfld{ifield});
76            if isequal(mncfld{ifield},'FW')
77                FocnOld = DataCorrections(FocnOld,mncfld{ifield},...
78                                          gravity,RAC,HFacC_Top);
79            end
80          Focn = (gWeight*FocnOld + FocnNew)/(gWeight+1);          Focn = (gWeight*FocnOld + FocnNew)/(gWeight+1);
81      end      end
82      fid=fopen([outDir,'/',FocnRoot,'.',fields{ifield},...      fid=fopen([outDir,'/',FocnRoot,'.',fields{ifield},...
# Line 78  for ifield = 1:length(mncfld) Line 85  for ifield = 1:length(mncfld)
85      fclose(fid);      fclose(fid);
86      dataFocn.(mncfld{ifield}) = Focn;      dataFocn.(mncfld{ifield}) = Focn;
87  end  end
88    % merccube_mod(XG,YG,dataFocn.HF);
89    % aaa = dataFocn.HF(:,:,end);
90    % sum(aaa(:).*RAC(:))./sum(RAC(:))
91    
92    
93  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94  %                         Prepare focn as needed                          %  %                         Prepare focn as needed                          %
95  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96    
97    % Must account for land!!!!
98  if ~FocnOnly  if ~FocnOnly
99            
100      years = [iitC:ditC:fitC].*dtC./31104000;      years = [iitC:ditC:fitC].*dtC./31104000;
101      yrfrc = rem(years,floor(years));      yrfrc = rem(years,floor(years));
102      yrfrc(yrfrc == 0) = 1;      yrfrc(yrfrc == 0) = 1;
103      uyrfrc = sort(unique(yrfrc));      uyrfrc = sort(unique(yrfrc));
104      length(yrfrc)      disp(['Cpl. data instances:  ',num2str(length(yrfrc))]);
105      length(uyrfrc)      disp(['Cpl. data per year forcing:  ',length(uyrfrc)]);
106            
107      for ifield = 1:length(mncfld)      for ifield = 1:length(mncfld)
108                    
109          % Make yearly forcing file.          % Make yearly forcing file (usually twelve months).
110          focnFull = dataCplSelect.([mncfld{ifield},'tave']);          focnFull = dataCplSelect.([mncfld{ifield},'tave']);
111          focn = NaN.*zeros([size(RAC,1),size(RAC,2),length(uyrfrc)]);          focn = NaN.*zeros([size(RAC,1),size(RAC,2),length(uyrfrc)]);
112          for itime = 1:length(uyrfrc)          for itime = 1:length(uyrfrc)
# Line 106  if ~FocnOnly Line 117  if ~FocnOnly
117              error('NaNs in ocean forcing files!!');              error('NaNs in ocean forcing files!!');
118          end          end
119                    
120          % Remove mean, add Focn          % 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          SubMean = mean(focn,3);          SubMean = mean(focn,3);
127          AddMean = dataFocn.(mncfld{ifield});          AddMean = dataFocn.(mncfld{ifield});
128          for itime = 1:size(focn,3)          for itime = 1:size(focn,3)
# Line 120  if ~FocnOnly Line 136  if ~FocnOnly
136              focn = focn.*iceMask;              focn = focn.*iceMask;
137          end          end
138                    
           
139          focnmean = mean(focn,3);          focnmean = mean(focn,3);
140          focnmean = sum(focnmean(:).*RAC(:))./sum(RAC(:));          focnmean = sum(focnmean(:).*RAC(:))./sum(RAC(:));
141          disp([mncfld{ifield},' mean:  ',num2str(focnmean)]);          disp([mncfld{ifield},' mean:  ',num2str(focnmean)]);
# Line 128  if ~FocnOnly Line 143  if ~FocnOnly
143          fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'w','b');          fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'w','b');
144          fwrite(fid,focn,'real*8');          fwrite(fid,focn,'real*8');
145          fclose(fid);          fclose(fid);
146            
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      end      end
151  end  end
152    
153    
154  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155  %                            Inspect Focn                                 %  %                            Inspect Focn                                 %
156  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 141  if InspectFocn Line 161  if InspectFocn
161                     '.bin.',num2str(period)],'r','b');                     '.bin.',num2str(period)],'r','b');
162          Focn = reshape(fread(fid,'real*8'),size(RAC));          Focn = reshape(fread(fid,'real*8'),size(RAC));
163          fclose(fid);          fclose(fid);
164          figure(ifield);          figure;
165          merccube(XG,YG,Focn); colorbar;          merccube_mod(XG,YG,Focn); colorbar;
166          title(mncfld{ifield});          title(mncfld{ifield});
167      end      end
168  end  end
169    
170  if Inspectfocn  if Inspectfocn
171      for ifield = 1:length(mncfld)      for ifield = 1:length(mncfld)
172          figure(ifield); orient tall;          figure; orient tall;
173          fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'r','b');          fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'r','b');
174          focn = reshape(fread(fid,'real*8'),[size(RAC),12]);          focn = reshape(fread(fid,'real*8'),[size(RAC),12]);
175          fclose(fid);          fclose(fid);
176          crange = [min(focn(:)),max(focn(:))];          crange = [min(focn(:)),max(focn(:))];
177          for ii = 1:size(focn,3)          for ii = 1:size(focn,3)
178              subplot(4,3,ii); merccube(XG,YG,focn(:,:,ii));              subplot(4,3,ii); merccube_mod(XG,YG,focn(:,:,ii));
179              shading flat; caxis(crange); colorbar;              shading flat; caxis(crange); colorbar;
180              title([mncfld{ifield},' month ',num2str(ii)]);              title([mncfld{ifield},' month ',num2str(ii)]);
181          end          end
# Line 165  end Line 185  end
185    
186  % Field dependent data corrections.  % Field dependent data corrections.
187  function CorrectedData = ...  function CorrectedData = ...
188    DataCorrections(DataToCorrect,CorrectionField,gravity,RAC)    DataCorrections(DataToCorrect,CorrectionField,gravity,RAC,HFacC_Top)
189        index = find(HFacC_Top);
190        correctdim = size(DataToCorrect);
191        DataToCorrect = DataToCorrect(:);
192      if     isequal(CorrectionField,'SIC')      if     isequal(CorrectionField,'SIC')
193          CorrectedData = gravity.*DataToCorrect;          CorrectedData(index) = gravity.*DataToCorrect(index);
194      elseif isequal(CorrectionField,'FW')      elseif isequal(CorrectionField,'FW')
195          offset = sum(DataToCorrect(:).*RAC(:))./sum(RAC(:));          DataToCorrect(index);
196          CorrectedData = DataToCorrect - offset;          offset = sum(DataToCorrect(index).*RAC(index))./sum(RAC(index));
197            CorrectedData(index) = DataToCorrect(index) - offset;
198      else      else
199          CorrectedData = DataToCorrect;          CorrectedData(index) = DataToCorrect(index);
200      end      end
201        CorrectedData = reshape(CorrectedData,correctdim);
202                    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22