--- MITgcm_contrib/enderton/PeriodicCoupling/CalcOcnForcingFields.m	2005/10/18 16:07:15	1.3
+++ MITgcm_contrib/enderton/PeriodicCoupling/CalcOcnForcingFields.m	2005/11/25 18:28:01	1.4
@@ -15,14 +15,15 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 datacpl = rdmnc(cplFiles,'iter','HFtave','TXtave',...
-                         'TYtave','FWtave','SICtave');
-datagrd = rdmnc(gridFiles,'rA','XG','YG','XC','YC');
+                'TYtave','FWtave','SICtave','T');
+datagrd = rdmnc(gridFiles,'rA','XG','YG','XC','YC','HFacC');
 RAC = datagrd.rA;
+HFacC = datagrd.HFacC;
 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');
-
+% save('DataPerCpl.mat','datacpl','RAC','XG','YG','XC','YC','HFacC');
+% load('DataPerCpl.mat');
+HFacC_Top = HFacC(:,:,1);
 
 % If a new focn is to be calculated, the old Focn must be loaded.
 if ~FocnOnly
@@ -35,6 +36,9 @@
         fclose(fid);
     end
 end
+% merccube_mod(XG,YG,dataFocnOld.HF);
+% aaa = dataFocnOld.HF;
+% sum(aaa(:).*RAC(:))./sum(RAC(:))
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -50,6 +54,9 @@
     dataCplSelect.([mncfld{ifield},'tave']) = ...
         datacpl.([mncfld{ifield},'tave'])(:,:,index);
 end
+% merccube_mod(XG,YG,dataCplSelect.HFtave(:,:,end));
+% aaa = dataCplSelect.HFtave(:,:,end);
+% sum(aaa(:).*RAC(:))./sum(RAC(:))
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -60,14 +67,15 @@
 % treatment of finding Focn.  Should compute monthly mean before.  Test
 % with various sorts of coupled field outputs.
 for ifield = 1:length(mncfld)
-    FocnNew = mean(dataCplSelect.([mncfld{ifield},'tave']),3);
-    FocnNew = DataCorrections(FocnNew,mncfld{ifield},gravity,RAC);
+    FocnNew = HFacC_Top.*mean(dataCplSelect.([mncfld{ifield},'tave']),3);
+    FocnNew = DataCorrections(FocnNew,mncfld{ifield},gravity,RAC,HFacC_Top);
     if FocnOnly
         Focn = FocnNew;
     else
-        FocnOld = dataFocnOld.(mncfld{ifield});
+        FocnOld = HFacC_Top.*dataFocnOld.(mncfld{ifield});
         if isequal(mncfld{ifield},'FW')
-            FocnOld = DataCorrections(FocnOld,mncfld{ifield},gravity,RAC);
+            FocnOld = DataCorrections(FocnOld,mncfld{ifield},...
+                                      gravity,RAC,HFacC_Top);
         end
         Focn = (gWeight*FocnOld + FocnNew)/(gWeight+1);
     end
@@ -77,24 +85,28 @@
     fclose(fid);
     dataFocn.(mncfld{ifield}) = Focn;
 end
+% merccube_mod(XG,YG,dataFocn.HF);
+% aaa = dataFocn.HF(:,:,end);
+% sum(aaa(:).*RAC(:))./sum(RAC(:))
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %                         Prepare focn as needed                          %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+% Must account for land!!!!
 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)
+    disp(['Cpl. data instances:  ',num2str(length(yrfrc))]);
+    disp(['Cpl. data per year forcing:  ',length(uyrfrc)]);
     
     for ifield = 1:length(mncfld)
         
-        % Make yearly forcing file.
+        % Make yearly forcing file (usually twelve months).
         focnFull = dataCplSelect.([mncfld{ifield},'tave']);
         focn = NaN.*zeros([size(RAC,1),size(RAC,2),length(uyrfrc)]);
         for itime = 1:length(uyrfrc)
@@ -105,7 +117,12 @@
             error('NaNs in ocean forcing files!!');
         end
         
-        % Remove mean, add Focn
+        % Account for bathymetry.
+        for itime = 1:size(focn,3)
+            focn(:,:,itime) = HFacC_Top.*focn(:,:,itime);
+        end
+        
+        % Remove time mean from focn, add Focn
         SubMean = mean(focn,3);
         AddMean = dataFocn.(mncfld{ifield});
         for itime = 1:size(focn,3)
@@ -119,7 +136,6 @@
             focn = focn.*iceMask;
         end
         
-        
         focnmean = mean(focn,3);
         focnmean = sum(focnmean(:).*RAC(:))./sum(RAC(:));
         disp([mncfld{ifield},' mean:  ',num2str(focnmean)]);
@@ -127,9 +143,14 @@
         fid=fopen([outDir,'/',focnRoot,'.',fields{ifield},'.bin'],'w','b');
         fwrite(fid,focn,'real*8');
         fclose(fid);
+        
+        aaa = mean(focn,3);
+        disp(['Number of land cells that do not have zero values:  ',...
+              num2str(length(find(aaa(find(~HFacC_Top))~=0)))]);
     end
 end
 
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %                            Inspect Focn                                 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -140,21 +161,21 @@
                    '.bin.',num2str(period)],'r','b');
         Focn = reshape(fread(fid,'real*8'),size(RAC));
         fclose(fid);
-        figure(ifield);
-        merccube(XG,YG,Focn); colorbar;
+        figure;
+        merccube_mod(XG,YG,Focn); colorbar;
         title(mncfld{ifield});
     end
 end
 
 if Inspectfocn
     for ifield = 1:length(mncfld)
-        figure(ifield); orient tall;
+        figure; 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));
+            subplot(4,3,ii); merccube_mod(XG,YG,focn(:,:,ii));
             shading flat; caxis(crange); colorbar;
             title([mncfld{ifield},' month ',num2str(ii)]);
         end
@@ -164,13 +185,18 @@
 
 % Field dependent data corrections.
 function CorrectedData = ...
-  DataCorrections(DataToCorrect,CorrectionField,gravity,RAC)
+  DataCorrections(DataToCorrect,CorrectionField,gravity,RAC,HFacC_Top)
+    index = find(HFacC_Top);
+    correctdim = size(DataToCorrect);
+    DataToCorrect = DataToCorrect(:);
     if     isequal(CorrectionField,'SIC')
-        CorrectedData = gravity.*DataToCorrect;
+        CorrectedData(index) = gravity.*DataToCorrect(index);
     elseif isequal(CorrectionField,'FW')
-        offset = sum(DataToCorrect(:).*RAC(:))./sum(RAC(:));
-        CorrectedData = DataToCorrect - offset;
+        DataToCorrect(index);
+        offset = sum(DataToCorrect(index).*RAC(index))./sum(RAC(index));
+        CorrectedData(index) = DataToCorrect(index) - offset;
     else
-        CorrectedData = DataToCorrect;
+        CorrectedData(index) = DataToCorrect(index);
     end
+    CorrectedData = reshape(CorrectedData,correctdim);