matlab_class/gcmfaces_diags/diags_budg_geothermal.m

function []=diags_budg_geothermal(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:
fileStartGeothFlux='geothermalFlux';
if strcmp(fileStart,'budg2d_snap_set2');
k=1; test3d=0;
elseif ~isempty(strfind(fileStart,'budg2d_snap_set3'));
k=str2num(fileStart(18:end)); test3d=0;
fileStartGeothFlux=['geothermalFlux_' fileStart(18:end)];
else;
k=1; test3d=1;
end;

%assemble the 3d geothermal flux field (flux entering a given layer at the underlying interface)
fld2d=read_bin([dirSnap 'geothermalFlux.bin']);
fld3d=0*mygrid.mskC;
fld3d=cat(3,fld3d,NaN*fld3d(:,:,1));
nr=length(mygrid.RC);
for kk=1:nr;
  tmp1=fld3d(:,:,kk);
  tmp2=fld3d(:,:,kk+1);
  tmp3=find(isnan(tmp2)&~isnan(tmp1)); 
  tmp1(tmp3)=fld2d(tmp3);
  tmp1(isnan(tmp1))=0;
  fld3d(:,:,kk)=tmp1;
end;

if test3d;
%output 3D field
fld2=fld3d;%maybe we need to switch by 1 level?
error('has not been tested yet');
%
else;
%integrate over underlying interfaces:
fld2=sum(fld3d(:,:,k:nr),3);
end;

%4) write to file:
fileMetaOld=[dirSnap fileList(ii+1).name(1:end-5) '.meta'];
fileIter=fileList(ii+1).name(length(fileStart)+2:1:end-5);
fileMetaNew=[dirMat 'BUDG/' fileStartGeothFlux '.' fileIter '.meta'];
%
fidMetaOld=fopen(fileMetaOld,'r');
fidMetaNew=fopen(fileMetaNew,'w');
test0=1;
while test0;
tmp1=fgetl(fidMetaOld);
test1=isempty(strfind(tmp1,'fldList'));
test2=isempty(strfind(tmp1,'nrecords'));
test3=isempty(strfind(tmp1,'nFlds'));
if test1&test2&test3;
  fprintf(fidMetaNew,[tmp1 '\n']);
elseif ~test3;
  fprintf(fidMetaNew,' nFlds = [    1 ];\n');
elseif ~test2;
  fprintf(fidMetaNew,' nrecords = [     1 ];\n');
else;
  test0=0;
end;
end;
fprintf(fidMetaNew,' fldList = {\n');
fprintf(fidMetaNew,'''geothFlux''\n');
fprintf(fidMetaNew,' };\n');
fclose(fidMetaOld);
fclose(fidMetaNew);
%
fileDataNew=[dirMat 'BUDG/' fileStartGeothFlux '.' fileIter '.data'];
write2file(fileDataNew,convert2gcmfaces(fld2),dataprec);
end;

1	function []=diags_budg_geothermal(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	fileStartGeothFlux='geothermalFlux';
25	if strcmp(fileStart,'budg2d_snap_set2');
26	k=1; test3d=0;
27	elseif ~isempty(strfind(fileStart,'budg2d_snap_set3'));
28	k=str2num(fileStart(18:end)); test3d=0;
29	fileStartGeothFlux=['geothermalFlux_' fileStart(18:end)];
30	else;
31	k=1; test3d=1;
32	end;
33
34	%assemble the 3d geothermal flux field (flux entering a given layer at the underlying interface)
35	fld2d=read_bin([dirSnap 'geothermalFlux.bin']);
36	fld3d=0*mygrid.mskC;
37	fld3d=cat(3,fld3d,NaN*fld3d(:,:,1));
38	nr=length(mygrid.RC);
39	for kk=1:nr;
40	tmp1=fld3d(:,:,kk);
41	tmp2=fld3d(:,:,kk+1);
42	tmp3=find(isnan(tmp2)&~isnan(tmp1));
43	tmp1(tmp3)=fld2d(tmp3);
44	tmp1(isnan(tmp1))=0;
45	fld3d(:,:,kk)=tmp1;
46	end;
47
48	if test3d;
49	%output 3D field
50	fld2=fld3d;%maybe we need to switch by 1 level?
51	error('has not been tested yet');
52	%
53	else;
54	%integrate over underlying interfaces:
55	fld2=sum(fld3d(:,:,k:nr),3);
56	end;
57
58	%4) write to file:
59	fileMetaOld=[dirSnap fileList(ii+1).name(1:end-5) '.meta'];
60	fileIter=fileList(ii+1).name(length(fileStart)+2:1:end-5);
61	fileMetaNew=[dirMat 'BUDG/' fileStartGeothFlux '.' fileIter '.meta'];
62	%
63	fidMetaOld=fopen(fileMetaOld,'r');
64	fidMetaNew=fopen(fileMetaNew,'w');
65	test0=1;
66	while test0;
67	tmp1=fgetl(fidMetaOld);
68	test1=isempty(strfind(tmp1,'fldList'));
69	test2=isempty(strfind(tmp1,'nrecords'));
70	test3=isempty(strfind(tmp1,'nFlds'));
71	if test1&test2&test3;
72	fprintf(fidMetaNew,[tmp1 '\n']);
73	elseif ~test3;
74	fprintf(fidMetaNew,' nFlds = [ 1 ];\n');
75	elseif ~test2;
76	fprintf(fidMetaNew,' nrecords = [ 1 ];\n');
77	else;
78	test0=0;
79	end;
80	end;
81	fprintf(fidMetaNew,' fldList = {\n');
82	fprintf(fidMetaNew,'''geothFlux''\n');
83	fprintf(fidMetaNew,' };\n');
84	fclose(fidMetaOld);
85	fclose(fidMetaNew);
86	%
87	fileDataNew=[dirMat 'BUDG/' fileStartGeothFlux '.' fileIter '.data'];
88	write2file(fileDataNew,convert2gcmfaces(fld2),dataprec);
89	end;
90