enderton/PeriodicCoupling/MakePcForcing.m

function MakePcForcing(cplFiles,atmFiles,gridFiles,prevFocnFiles,...
                       period,outDir,FocnRoot,focnRoot,iterO,iterA,...
                       dtC,dtA,gWeight,FocnOnly,InspectFocn,Inspectfocn)

% Make fields for forcing ocean component of periodic coupling.

rhofw = 1000.;

fieldscpl = {'TXtave'  ,'TYtave'  ,'FWtave' ,'HFtave'  };
fieldsatm = {'THETA'   ,'ALBVISDF','RADSWG'  ,...
             'DWNLWG'  ,'TS'      ,'QS'      ,...
             'WINDS'   ,'PRECON'  ,'PRECLS'};
fieldsall = {'TXtave'  ,'TYtave'  ,'FWtave'  ,'HFtave'  ,...
             'THETA'   ,'DWNLWG'  ,'TS'      ,'QS'      ,...
             'WINDS'   ,'DWNSWG'  ,'PREC'  };

warning off MATLAB:MKDIR:DirectoryExists;
mkdir(outDir);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 Load Data                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

datacpl = rdmnc(cplFiles,'iter','HFtave','TXtave','TYtave','FWtave','T');
dataatm = rdmnc(atmFiles,'iter','THETA'   ,'ALBVISDF','RADSWG'  ,...
                                'DWNLWG'  ,'TS'      ,'QS'      ,...
                                'WINDS'   ,'PRECON'  ,'PRECLS','T');
datagrd = rdmnc(gridFiles,'rA','XG','YG','XC','YC','HFacC');
% save('DataPerCpl.mat','datacpl','dataatm','datagrd');
% load('DataPerCpl.mat');
nc = size(datagrd.rA,2);
RAC = datagrd.rA;
mask = datagrd.HFacC(:,:,1);
imask = find(mask);
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));

% If a new focn is to be calculated, the old Focn must be loaded.
if ~FocnOnly
    for ifield=1:length(fieldsall)
        fid=fopen([prevFocnFiles,'/',FocnRoot,...
                   '.',fieldsall{ifield},'.bin.',...
                   num2str(period-1)],'r','b');
        dataFocnOld.(fieldsall{ifield}) = ...
            reshape(fread(fid,'real*8'),size(RAC));
        fclose(fid);
    end
end
% field = 'HFtave';
% merccube_mod(XG,YG,dataFocnOld.(field));
% aaa = dataFocnOld.(field)(:,:,end);
% sum(aaa(:).*RAC(:))./sum(RAC(:))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      Select Coupled Field Data                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isequal(datacpl.iter,dataatm.iter./(dtC./dtA))
    disp('datacpl.iter');           datacpl.iter
    disp('dataatm.iter/dt_ratio');  dataatm.iter./(dtC./dtA)
    error('Iterations for coupled and atmospheric data do not match.');
end
test   = ismember(iterO,datacpl.iter);
indexO = ismember(datacpl.iter,iterO);
if ~isempty(find(test==0))
    error('Specified indecies not present in cpl_tave.* files.');
else
    for ifield = 1:length(fieldscpl)
        data.(fieldscpl{ifield}) = ...
            datacpl.(fieldscpl{ifield})(1:6*nc,1:nc,indexO);
    end
end
test   = ismember(iterA,dataatm.iter);
indexA = ismember(dataatm.iter,iterA);
if ~isempty(find(test==0))
    error('Specified indecies not present in DiagAtmForcing.* files.');
else
    for ifield = 1:length(fieldsatm)
        data.(fieldsatm{ifield}) = ...
            dataatm.(fieldsatm{ifield})(1:6*nc,1:nc,indexA);
    end
end
% field = 'PREC';
% merccube_mod(XG,YG,dataCplSelect.(field));
% aaa = dataCplSelect.(field)(:,:,end);
% sum(aaa(:).*RAC(:))./sum(RAC(:))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Prepare Focn                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% What if it is not in year increments?  Need a more sophisticated
% treatment of finding Focn.  Should compute monthly mean before.  Test
% with various sorts of coupled field outputs.
for ifield = 1:length(fieldsall)
    if isequal(fieldsall{ifield},'PREC')
        FocnNew = mask.*( mean(data.PRECON,3) ...
                         +mean(data.PRECLS,3))./(1000.*rhofw);
    elseif isequal(fieldsall{ifield},'DWNSWG')
        FocnNew = mask.*mean(data.RADSWG,3)./(1-mean(data.ALBVISDF,3));
    else
        FocnNew = mask.*mean(data.(fieldsall{ifield}),3);
    end
    if FocnOnly
        Focn = FocnNew;
    else
        FocnOld = mask.*dataFocnOld.(fieldsall{ifield});
        Focn = (gWeight*FocnOld + FocnNew)/(gWeight+1);
    end
    if isequal(fieldsall{ifield},'FWtave')
        offset = sum(Focn(imask).*RAC(imask))./sum(RAC(imask));
        Focn(imask) = Focn(imask) - offset;
    end
    dataFocn.(fieldsall{ifield}) = Focn;
    fid=fopen([outDir,'/',FocnRoot,'.',...
               fieldsall{ifield},'.bin.',...
               num2str(period)],'w','b');
    fwrite(fid,Focn,'real*8');
    fclose(fid);
end
% field = 'PREC';
% merccube_mod(XG,YG,dataFocn.(field));
% aaa = dataFocn.(field)(:,:,end);
% sum(aaa(:).*RAC(:))./sum(RAC(:))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         Prepare focn as needed                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~FocnOnly
    
    years = iterO.*dtC./31104000;
    yrfrc = rem(years,floor(years));
    yrfrc(yrfrc == 0) = 1;
    uyrfrc = sort(unique(yrfrc));
    disp(['Cpl. data instances:  ',num2str(length(yrfrc))]);
    disp(['Cpl. data per year forcing:  ',num2str(length(uyrfrc))]);
    
    for ifield = 1:length(fieldsall)
        
        % Make yearly forcing file (usually twelve months).
        if isequal(fieldsall{ifield},'PREC')
            focnFull = (data.PRECON+data.PRECLS)./(1000.*rhofw);
        elseif isequal(fieldsall{ifield},'DWNSWG')
            focnFull = (data.RADSWG)./(1-data.ALBVISDF);
        else
            focnFull = data.(fieldsall{ifield});
        end
        focn = NaN.*zeros([size(RAC,1),size(RAC,2),length(uyrfrc)]);
        for itime = 1:length(uyrfrc)
            index = find(uyrfrc(itime) == yrfrc);
            focn(:,:,itime) = mean(focnFull(:,:,index),3);
        end
        if ~isempty(find(isnan(focn(:))))
            error('NaNs in ocean forcing files!!');
        end
        
        % Account for bathymetry.
        for itime = 1:size(focn,3)
            focn(:,:,itime) = mask.*focn(:,:,itime);
        end
        
        % Remove time mean from focn, add Focn
        SubMean = mean(focn,3);
        AddMean = dataFocn.(fieldsall{ifield});
        for itime = 1:size(focn,3)
            focn(:,:,itime) = focn(:,:,itime) - SubMean + AddMean;
        end
        
        focnmean = mean(focn,3);
        focnmean = sum(focnmean(:).*RAC(:))./sum(RAC(:));
        disp([fieldsall{ifield},' mean:  ',num2str(focnmean)]);
        
        fid=fopen([outDir,'/',focnRoot,'.',...
                   fieldsall{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(~mask))~=0)))]);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Inspect Focn                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if InspectFocn
    for ifield = 1:length(fieldsall)
        fid=fopen([outDir,'/',FocnRoot,'.',...
                   fieldsall{ifield},'.bin.',...
                   num2str(period)],'r','b');
        Focn = reshape(fread(fid,'real*8'),size(RAC));
        fclose(fid);
        figure; merccube_mod(XG,YG,Focn); colorbar;
        title(fieldsall{ifield});
    end
end

if Inspectfocn
    for ifield = 1:length(fieldsall)
        fid=fopen([outDir,'/',focnRoot,'.',...
                   fieldsall{ifield},'.bin'],...
                   'r','b');
        focn = reshape(fread(fid,'real*8'),[size(RAC),12]);
        fclose(fid);
        figure;
        crange = [min(focn(:)),max(focn(:))];
        for ii = 1:size(focn,3)
            subplot(4,3,ii); merccube_mod(XG,YG,focn(:,:,ii));
            shading flat; caxis(crange); colorbar;
            title([fieldsall{ifield},' month ',num2str(ii)]);
        end
    end
end
1	enderton	1.2	function MakePcForcing(cplFiles,atmFiles,gridFiles,prevFocnFiles,...
2	enderton	1.3	period,outDir,FocnRoot,focnRoot,iterO,iterA,...
3			dtC,dtA,gWeight,FocnOnly,InspectFocn,Inspectfocn)
4	enderton	1.1
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	enderton	1.4	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	enderton	1.1	% save('DataPerCpl.mat','datacpl','dataatm','datagrd');
31	enderton	1.4	% load('DataPerCpl.mat');
32	enderton	1.1	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	enderton	1.3	if ~isequal(datacpl.iter,dataatm.iter./(dtC./dtA))
61			disp('datacpl.iter'); datacpl.iter
62			disp('dataatm.iter/dt_ratio'); dataatm.iter./(dtC./dtA)
63	enderton	1.2	error('Iterations for coupled and atmospheric data do not match.');
64			end
65	enderton	1.3	test = ismember(iterO,datacpl.iter);
66			indexO = ismember(datacpl.iter,iterO);
67	enderton	1.1	if ~isempty(find(test==0))
68			error('Specified indecies not present in cpl_tave.* files.');
69	enderton	1.3	else
70			for ifield = 1:length(fieldscpl)
71			data.(fieldscpl{ifield}) = ...
72			datacpl.(fieldscpl{ifield})(1:6*nc,1:nc,indexO);
73			end
74	enderton	1.1	end
75	enderton	1.3	test = ismember(iterA,dataatm.iter);
76			indexA = ismember(dataatm.iter,iterA);
77			if ~isempty(find(test==0))
78			error('Specified indecies not present in DiagAtmForcing.* files.');
79			else
80			for ifield = 1:length(fieldsatm)
81			data.(fieldsatm{ifield}) = ...
82			dataatm.(fieldsatm{ifield})(1:6*nc,1:nc,indexA);
83			end
84	enderton	1.1	end
85			% field = 'PREC';
86			% merccube_mod(XG,YG,dataCplSelect.(field));
87			% aaa = dataCplSelect.(field)(:,:,end);
88			% sum(aaa(:).*RAC(:))./sum(RAC(:))
89
90
91			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
92			% Prepare Focn %
93			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94
95			% What if it is not in year increments? Need a more sophisticated
96			% treatment of finding Focn. Should compute monthly mean before. Test
97			% with various sorts of coupled field outputs.
98			for ifield = 1:length(fieldsall)
99			if isequal(fieldsall{ifield},'PREC')
100			FocnNew = mask.*( mean(data.PRECON,3) ...
101			+mean(data.PRECLS,3))./(1000.*rhofw);
102			elseif isequal(fieldsall{ifield},'DWNSWG')
103			FocnNew = mask.*mean(data.RADSWG,3)./(1-mean(data.ALBVISDF,3));
104			else
105			FocnNew = mask.*mean(data.(fieldsall{ifield}),3);
106			end
107			if FocnOnly
108			Focn = FocnNew;
109			else
110			FocnOld = mask.*dataFocnOld.(fieldsall{ifield});
111			Focn = (gWeight*FocnOld + FocnNew)/(gWeight+1);
112			end
113			if isequal(fieldsall{ifield},'FWtave')
114			offset = sum(Focn(imask).*RAC(imask))./sum(RAC(imask));
115			Focn(imask) = Focn(imask) - offset;
116			end
117			dataFocn.(fieldsall{ifield}) = Focn;
118			fid=fopen([outDir,'/',FocnRoot,'.',...
119			fieldsall{ifield},'.bin.',...
120			num2str(period)],'w','b');
121			fwrite(fid,Focn,'real*8');
122			fclose(fid);
123			end
124			% field = 'PREC';
125			% merccube_mod(XG,YG,dataFocn.(field));
126			% aaa = dataFocn.(field)(:,:,end);
127			% sum(aaa(:).*RAC(:))./sum(RAC(:))
128
129
130			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131			% Prepare focn as needed %
132			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
133
134			if ~FocnOnly
135
136	enderton	1.3	years = iterO.*dtC./31104000;
137	enderton	1.1	yrfrc = rem(years,floor(years));
138			yrfrc(yrfrc == 0) = 1;
139			uyrfrc = sort(unique(yrfrc));
140			disp(['Cpl. data instances: ',num2str(length(yrfrc))]);
141			disp(['Cpl. data per year forcing: ',num2str(length(uyrfrc))]);
142
143			for ifield = 1:length(fieldsall)
144
145			% Make yearly forcing file (usually twelve months).
146			if isequal(fieldsall{ifield},'PREC')
147			focnFull = (data.PRECON+data.PRECLS)./(1000.*rhofw);
148			elseif isequal(fieldsall{ifield},'DWNSWG')
149			focnFull = (data.RADSWG)./(1-data.ALBVISDF);
150			else
151			focnFull = data.(fieldsall{ifield});
152			end
153			focn = NaN.*zeros([size(RAC,1),size(RAC,2),length(uyrfrc)]);
154			for itime = 1:length(uyrfrc)
155			index = find(uyrfrc(itime) == yrfrc);
156			focn(:,:,itime) = mean(focnFull(:,:,index),3);
157			end
158			if ~isempty(find(isnan(focn(:))))
159			error('NaNs in ocean forcing files!!');
160			end
161
162			% Account for bathymetry.
163			for itime = 1:size(focn,3)
164			focn(:,:,itime) = mask.*focn(:,:,itime);
165			end
166
167			% Remove time mean from focn, add Focn
168			SubMean = mean(focn,3);
169			AddMean = dataFocn.(fieldsall{ifield});
170			for itime = 1:size(focn,3)
171			focn(:,:,itime) = focn(:,:,itime) - SubMean + AddMean;
172			end
173
174			focnmean = mean(focn,3);
175			focnmean = sum(focnmean(:).*RAC(:))./sum(RAC(:));
176			disp([fieldsall{ifield},' mean: ',num2str(focnmean)]);
177
178			fid=fopen([outDir,'/',focnRoot,'.',...
179			fieldsall{ifield},'.bin'],'w','b');
180			fwrite(fid,focn,'real*8');
181			fclose(fid);
182
183			aaa = mean(focn,3);
184			disp(['Number of land cells that do not have zero values: ',...
185			num2str(length(find(aaa(find(~mask))~=0)))]);
186			end
187			end
188
189
190			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
191			% Inspect Focn %
192			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
193
194			if InspectFocn
195			for ifield = 1:length(fieldsall)
196			fid=fopen([outDir,'/',FocnRoot,'.',...
197			fieldsall{ifield},'.bin.',...
198			num2str(period)],'r','b');
199			Focn = reshape(fread(fid,'real*8'),size(RAC));
200			fclose(fid);
201			figure; merccube_mod(XG,YG,Focn); colorbar;
202			title(fieldsall{ifield});
203			end
204			end
205
206			if Inspectfocn
207			for ifield = 1:length(fieldsall)
208			fid=fopen([outDir,'/',focnRoot,'.',...
209			fieldsall{ifield},'.bin'],...
210			'r','b');
211			focn = reshape(fread(fid,'real*8'),[size(RAC),12]);
212			fclose(fid);
213			figure;
214			crange = [min(focn(:)),max(focn(:))];
215			for ii = 1:size(focn,3)
216			subplot(4,3,ii); merccube_mod(XG,YG,focn(:,:,ii));
217			shading flat; caxis(crange); colorbar;
218			title([fieldsall{ifield},' month ',num2str(ii)]);
219			end
220			end
221			end