matlab_class/gcmfaces_diags/diags_diff_snapshots.m

function []=diags_diff_snapshots(dirSnap,dirMat,fileStart);
%object: compute time derivatives between snapshots that 
%   will be compared with time mean flux terms in budgets
%input:     dirSnap is the directory containing the binary snapshots
%           fileStart is the root of file names (e.g. budg2d_snap_set1)
%result:    create e.g. rate_budg2d_snap_set1 files

gcmfaces_global; global myparms;

fileList=dir([dirSnap fileStart '.*.data']);
for ii=1:length(fileList)-1;

%1) get the time of fld0 & fld1, & the data precision
fileName0=fullfile(dirSnap,fileList(ii).name(1:end-5));
meta0=rdmds_meta(fileName0);
fileName1=fullfile(dirSnap,filesep,fileList(ii+1).name(1:end-5));
meta1=rdmds_meta(fileName1);
%
time1=meta1.timeInterval;
time0=meta0.timeInterval;
dataprec=str2num(meta1.dataprec(end-1:end));

%2) get the binary data:
fld0=rdmds(fileName0);
fld1=rdmds(fileName1);

test3d=(length(size(fld0))==4);
if test3d&(myparms.useNLFS==2);
%3D diagnostics are multiplied by DRF*hFac*ETAN
%(2D diagnostics are expectedly vertically integrated by MITgcm)
  for jj=0:1;
    %get etan
    etanName=['budg2d_snap_set1' fileList(ii+jj).name(end-15:end-5)];
    etanName=fullfile(dirSnap,filesep,etanName);
    etan=rdmds2gcmfaces(etanName,'rec',1);
    %compute time variable thickness
    tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
    tmp2=tmp1./mk3D(mygrid.Depth,tmp1);
    tmp2=tmp2.*mk3D(etan,tmp1);
    drf=convert2gcmfaces(tmp1+tmp2);
    %apply to scale fld
    eval(['fld=fld' num2str(jj) ';']);
    n4=size(fld,4);
    drf=repmat(drf,[1 1 1 n4]);
    fld=fld.*drf;
    eval(['fld' num2str(jj) '=fld;']); 
  end;
elseif test3d;
  %the non rstar case remains to be treated
  error('missing implementation of diags_diff_snapshots\n');
end;

%3) compute the tendency term:
fld2=(fld1-fld0)/(time1-time0);

%4) write to file:
fileMetaOld=[dirSnap fileList(ii+1).name(1:end-5) '.meta'];
fileMetaNew=[dirMat 'BUDG/rate_' fileList(ii+1).name(1:end-5) '.meta'];
copyfile(fileMetaOld , fileMetaNew);
fileDataNew=[dirMat 'BUDG/rate_' fileList(ii+1).name(1:end-5) '.data'];
write2file(fileDataNew,fld2,dataprec);

end;

1	function []=diags_diff_snapshots(dirSnap,dirMat,fileStart);
2	%object: compute time derivatives between snapshots that
3	% will be compared with time mean flux terms in budgets
4	%input: dirSnap is the directory containing the binary snapshots
5	% fileStart is the root of file names (e.g. budg2d_snap_set1)
6	%result: create e.g. rate_budg2d_snap_set1 files
7
8	gcmfaces_global; global myparms;
9
10	fileList=dir([dirSnap fileStart '.*.data']);
11	for ii=1:length(fileList)-1;
12
13	%1) get the time of fld0 & fld1, & the data precision
14	fileName0=fullfile(dirSnap,fileList(ii).name(1:end-5));
15	meta0=rdmds_meta(fileName0);
16	fileName1=fullfile(dirSnap,filesep,fileList(ii+1).name(1:end-5));
17	meta1=rdmds_meta(fileName1);
18	%
19	time1=meta1.timeInterval;
20	time0=meta0.timeInterval;
21	dataprec=str2num(meta1.dataprec(end-1:end));
22
23	%2) get the binary data:
24	fld0=rdmds(fileName0);
25	fld1=rdmds(fileName1);
26
27	test3d=(length(size(fld0))==4);
28	if test3d&(myparms.useNLFS==2);
29	%3D diagnostics are multiplied by DRFhFacETAN
30	%(2D diagnostics are expectedly vertically integrated by MITgcm)
31	for jj=0:1;
32	%get etan
33	etanName=['budg2d_snap_set1' fileList(ii+jj).name(end-15:end-5)];
34	etanName=fullfile(dirSnap,filesep,etanName);
35	etan=rdmds2gcmfaces(etanName,'rec',1);
36	%compute time variable thickness
37	tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
38	tmp2=tmp1./mk3D(mygrid.Depth,tmp1);
39	tmp2=tmp2.*mk3D(etan,tmp1);
40	drf=convert2gcmfaces(tmp1+tmp2);
41	%apply to scale fld
42	eval(['fld=fld' num2str(jj) ';']);
43	n4=size(fld,4);
44	drf=repmat(drf,[1 1 1 n4]);
45	fld=fld.*drf;
46	eval(['fld' num2str(jj) '=fld;']);
47	end;
48	elseif test3d;
49	%the non rstar case remains to be treated
50	error('missing implementation of diags_diff_snapshots\n');
51	end;
52
53	%3) compute the tendency term:
54	fld2=(fld1-fld0)/(time1-time0);
55
56	%4) write to file:
57	fileMetaOld=[dirSnap fileList(ii+1).name(1:end-5) '.meta'];
58	fileMetaNew=[dirMat 'BUDG/rate_' fileList(ii+1).name(1:end-5) '.meta'];
59	copyfile(fileMetaOld , fileMetaNew);
60	fileDataNew=[dirMat 'BUDG/rate_' fileList(ii+1).name(1:end-5) '.data'];
61	write2file(fileDataNew,fld2,dataprec);
62
63	end;
64