matlab_class/gcmfaces_diags/diags_diff_snapshots.m

function []=diags_diff_snapshots(dirModel,dirMat,fileStart);
%object: compute time derivatives between snapshots that 
%   will be compared with time mean flux terms in budgets
%input:     dirModel is the model directory
%           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([dirModel 'diags/BUDG/' fileStart '.*.data']);
for ii=1:length(fileList)-1;

%1) get the time of fld0 & fld1, & the data precision
fileName0=fullfile(dirModel,'diags',filesep,'BUDG',filesep,fileList(ii).name(1:end-5));
meta0=rdmds_meta(fileName0);
fileName1=fullfile(dirModel,'diags',filesep,'BUDG',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(dirModel,'diags',filesep,'BUDG',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=[dirModel 'diags/BUDG/' 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	gforget	1.1	function []=diags_diff_snapshots(dirModel,dirMat,fileStart);
2			%object: compute time derivatives between snapshots that
3			% will be compared with time mean flux terms in budgets
4			%input: dirModel is the model directory
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	gforget	1.3	gcmfaces_global; global myparms;
9
10	gforget	1.1	fileList=dir([dirModel 'diags/BUDG/' fileStart '.*.data']);
11			for ii=1:length(fileList)-1;
12
13			%1) get the time of fld0 & fld1, & the data precision
14	gforget	1.3	fileName0=fullfile(dirModel,'diags',filesep,'BUDG',filesep,fileList(ii).name(1:end-5));
15			meta0=rdmds_meta(fileName0);
16			fileName1=fullfile(dirModel,'diags',filesep,'BUDG',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	gforget	1.1
23			%2) get the binary data:
24	gforget	1.3	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(dirModel,'diags',filesep,'BUDG',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	gforget	1.1
53			%3) compute the tendency term:
54			fld2=(fld1-fld0)/(time1-time0);
55
56			%4) write to file:
57			fileMetaOld=[dirModel 'diags/BUDG/' fileList(ii+1).name(1:end-5) '.meta'];
58			fileMetaNew=[dirMat 'BUDG/rate_' fileList(ii+1).name(1:end-5) '.meta'];
59	gforget	1.2	copyfile(fileMetaOld , fileMetaNew);
60	gforget	1.1	fileDataNew=[dirMat 'BUDG/rate_' fileList(ii+1).name(1:end-5) '.data'];
61			write2file(fileDataNew,fld2,dataprec);
62
63			end;
64