matlab_class/gcmfaces_IO/prep2nctiles.m


choiceStruct=2;

%===============

if choiceStruct==1;
    structIn=[];
    %
    structIn.vars=mygrid;
    structIn.vars.RF=structIn.vars.RF(1:50);
    %
    structIn.descr={'C-grid parameters (see MITgcm documentation for details).'};
    %
    vars=[]; nv=length(vars)+1;
    vars(nv).fldName='XC'; vars(nv).longName='longitude'; vars(nv).units='degrees_east'; nv=length(vars)+1;
    vars(nv).fldName='YC'; vars(nv).longName='latitude'; vars(nv).units='degrees_north'; nv=length(vars)+1;
    vars(nv).fldName='XG'; vars(nv).longName='longitude'; vars(nv).units='degrees_east'; nv=length(vars)+1;
    vars(nv).fldName='YG'; vars(nv).longName='latitude'; vars(nv).units='degrees_north'; nv=length(vars)+1;
    vars(nv).fldName='RAC'; vars(nv).longName='grid cell area'; vars(nv).units='m^2'; nv=length(vars)+1;
    vars(nv).fldName='RAZ'; vars(nv).longName='grid cell area'; vars(nv).units='m^2'; nv=length(vars)+1;
    vars(nv).fldName='DXC'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='DYC'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='DXG'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='DYG'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='hFacC'; vars(nv).longName='fractional thickness'; vars(nv).units='1'; nv=length(vars)+1;
    vars(nv).fldName='hFacW'; vars(nv).longName='fractional thickness'; vars(nv).units='1'; nv=length(vars)+1;
    vars(nv).fldName='hFacS'; vars(nv).longName='fractional thickness'; vars(nv).units='1'; nv=length(vars)+1;
    vars(nv).fldName='Depth'; vars(nv).longName='sea floor depth'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='AngleCS'; vars(nv).longName='grid orientation (cosine)'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='AngleSN'; vars(nv).longName='grid orientation (sine)'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='RC'; vars(nv).longName='vertical coordinate'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='RF'; vars(nv).longName='vertical coordinate'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='DRC'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
    vars(nv).fldName='DRF'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
    %
    structIn.defs=vars;
    %
    struct2nctiles('release1/','GRID',structIn,[90 90]);
end;

%===============

if choiceStruct==2;

  listBudgs={'budgMo','budgHo','budgSo','budgMi','budgHi','budgSi'};

  %directories and snapshot
  dirIn='r4it11.c65i/';
  dirMat='mat_budg3d/';
  fldName='snapshot';
  timName='0000000732'; suff='initial';
% timName='0000174564'; suff='final';
  t0Name=num2str(3600*str2num(timName));

  %load myparms
  eval(['load ' dirIn dirMat 'diags_grid_parms.mat myparms;']);

  %read THETA, SALT
  diagName=['budg3d_snap_set1.' timName];
  diagName=fullfile(dirIn,'diags',filesep,'BUDG',filesep,diagName);
  THETA=rdmds2gcmfaces(diagName,'rec',1);
  SALT=rdmds2gcmfaces(diagName,'rec',2);
  ONE=1*(SALT>0);

  %compute cell thickness
  etanName=['budg2d_snap_set1.' timName];
  etanName=fullfile(dirIn,'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);
  tmp3=mk3D(mygrid.RAC,tmp1);
  vol=tmp3.*(tmp1+tmp2);

  %load ice and snow volume/m2
  SIheff=rdmds2gcmfaces(etanName,'rec',2);
  SIsnow=rdmds2gcmfaces(etanName,'rec',3);

  for ii=1:length(listBudgs);
    budgName=listBudgs{ii};
    [budgName ' -- ' fldName ' -- ' timName]

    %directories and snapshot
    dirOut=fullfile(dirIn,'nctiles_budg',filesep);
    if ~isdir(dirOut); mkdir(dirOut); end;
    dirOut=fullfile(dirOut,budgName,filesep);
    if ~isdir(dirOut); mkdir(dirOut); end;

    switch budgName;
    case 'budgMo'; fld=ONE; fac=myparms.rhoconst;
    case 'budgHo'; fld=THETA; fac=myparms.rcp;
    case 'budgSo'; fld=SALT; fac=myparms.rhoconst;
    case 'budgMi'; FACheff=myparms.rhoi; FACsnow=myparms.rhosn;
    case 'budgHi'; FACheff=-myparms.flami*myparms.rhoi; FACsnow=-myparms.flami*myparms.rhosn;
    case 'budgSi'; FACheff=myparms.SIsal0*myparms.rhoi; FACsnow=0;
    end;

    %compute extensive snapshot
    structIn=[];
    if budgName(end)=='o';
      structIn.vars.snapshot=fac*fld.*vol;
    else;
      structIn.vars.snapshot=mygrid.RAC.*(FACheff*SIheff+FACsnow*SIsnow);
    end;

    %general description
    specs=[];
    if strcmp(budgName(5),'H');
        structIn.descr={'Heat budget in extensive form (in Watt, on C-Grid)'};
        specs.units='J';
        specs.name='heat content';
    elseif strcmp(budgName(5),'M');
        structIn.descr={'Mass budget in extensive form (in kg/s, on C-Grid)'};
        specs.units='kg';
        specs.name='mass content';
    elseif strcmp(budgName(5),'S');
        structIn.descr={'Salt budget in extensive form (in g/s, on C-Grid)'};
        specs.units='g';
        specs.name='salt content';
    end;

    %variables description
    vars=[]; nv=length(vars)+1;
    vars(nv).fldName=fldName; vars(nv).units=specs.units;
    if strcmp(budgName(6),'o'); vars(nv).longName=['ocean ' specs.name ' snapshot'];
    elseif strcmp(budgName(6),'i'); vars(nv).longName=['ice+snow ' specs.name ' snapshot'];
    end;
    vars(nv).longName=[vars(nv).longName ' at ' suff ' time=' t0Name 's'];
    %
    structIn.defs=vars;

    %create file
    tic; struct2nctiles(dirIn,fldName,structIn,[90 90]); toc;
    eval(['!mv ' pwd filesep dirIn filesep 'tmp_nctiles' filesep fldName ' ' pwd filesep dirOut suff]);

  end;%for ii=1:length(listBudgs);

end;

%===============

if choiceStruct==3;

  %budgName='budgHo';

  if budgName(end)=='o'; listFlds={'tend','trU','trV','trWtop'};
  else; listFlds={'tend','trU','trV','trWtop','trWbot'};
  end;

  for ii=1:length(listFlds);
    fldName=listFlds{ii};
    %fldName='trWtop';
    [budgName ' -- ' fldName]

    %directories
    dirIn='r4it11.c65i/';
    dirMat='mat_budg3d/';
    dirOut=fullfile(dirIn,'nctiles_budg',filesep);
    if ~isdir(dirOut); mkdir(dirOut); end;
    dirOut=fullfile(dirOut,budgName,filesep);
    if ~isdir(dirOut); mkdir(dirOut); end;

    %load variable
    eval(['load ' dirIn dirMat 'diags_grid_parms.mat myparms;']);
    dirMat=[dirIn dirMat 'diags_set_' budgName '/'];
    fileMat=[budgName '_*.mat'];
    tic; structIn.vars=diags_read_from_mat(dirMat,fileMat,fldName); toc;

    %time vectors
    [listTimes]=diags_list_times({[dirIn 'diags/BUDG/']},{'budg2d_hflux_set1'});
    structIn.vars.t0=3600*listTimes(1:end-2);
    structIn.vars.t1=3600*listTimes(2:end-1);

    %rename trWtop as trW if adequate
    if strcmp(fldName,'trWtop')&(budgName(end)=='o');
%note: since geothermal heating was added we likely need to add a nr+1 level to trW
%  using trWbot(:,:,nr) or output trWtop and trWbot themselves (as done for seaice)
      structIn.vars=setfield(structIn.vars,'trW',structIn.vars.trWtop);
      structIn.vars=rmfield(structIn.vars,'trWtop');
      structIn.vars.listDiags={'trW'};
      fldName='trW';
    end;

    %switch back to upward convention
    if strcmp(fldName(1:3),'trW');
%note: the following was needed in preparing the 03-Feb-2015 nctiles_budget (3d) files 
%   but is no longer needed since the trW sign convention in gcmfaces_diags was 
%   reversed to positive upward for consistency with MITgcm convention.
      tmp1=getfield(structIn.vars,fldName);
      structIn.vars=setfield(structIn.vars,fldName,-tmp1);
    end;

    %general description
    tmp1=diags_read_from_mat(dirMat,fileMat,'specs',1);
    specs=tmp1.specs;
    if strcmp(specs.units,'W');
        structIn.descr={'Heat budget in extensive form (in Watt, on C-Grid)'};
    elseif strcmp(specs.units,'kg/s');
        structIn.descr={'Mass budget in extensive form (in kg/s, on C-Grid)'};
    elseif strcmp(specs.units,'g/s');
        structIn.descr={'Salt budget in extensive form (in g/s, on C-Grid)'};
    else;
        error('unknown budget');
    end;

    %variables description
    vars=[]; nv=length(vars)+1;
    vars(nv).fldName=fldName; vars(nv).units=specs.units;
    switch fldName;
    case 'tend'; vars(nv).longName='tendency term';
    case 'trU'; vars(nv).longName='horizontal transport (U)';
    case 'trV'; vars(nv).longName='horizontal transport (V)';
    case 'trW'; vars(nv).longName='upward vertical transport (W)';
    case 'trWtop'; vars(nv).longName='upward vertical transport (W)';
    case 'trWbot'; vars(nv).longName='upward vertical transport (W)';
    end;
    nv=length(vars)+1;
    vars(nv).fldName='t0'; vars(nv).longName='initial time'; vars(nv).units='s'; nv=length(vars)+1;
    vars(nv).fldName='t1'; vars(nv).longName='final time'; vars(nv).units='s'; nv=length(vars)+1;
    %
    structIn.defs=vars;

    %create file
    tic; struct2nctiles(dirIn,fldName,structIn,[90 90]); toc;
    eval(['!mv ' pwd filesep dirIn filesep 'tmp_nctiles' filesep fldName ' ' pwd filesep dirOut]);

  end;%for fldName=listFlds;

end;

1
2	choiceStruct=2;
3
4	%===============
5
6	if choiceStruct==1;
7	structIn=[];
8	%
9	structIn.vars=mygrid;
10	structIn.vars.RF=structIn.vars.RF(1:50);
11	%
12	structIn.descr={'C-grid parameters (see MITgcm documentation for details).'};
13	%
14	vars=[]; nv=length(vars)+1;
15	vars(nv).fldName='XC'; vars(nv).longName='longitude'; vars(nv).units='degrees_east'; nv=length(vars)+1;
16	vars(nv).fldName='YC'; vars(nv).longName='latitude'; vars(nv).units='degrees_north'; nv=length(vars)+1;
17	vars(nv).fldName='XG'; vars(nv).longName='longitude'; vars(nv).units='degrees_east'; nv=length(vars)+1;
18	vars(nv).fldName='YG'; vars(nv).longName='latitude'; vars(nv).units='degrees_north'; nv=length(vars)+1;
19	vars(nv).fldName='RAC'; vars(nv).longName='grid cell area'; vars(nv).units='m^2'; nv=length(vars)+1;
20	vars(nv).fldName='RAZ'; vars(nv).longName='grid cell area'; vars(nv).units='m^2'; nv=length(vars)+1;
21	vars(nv).fldName='DXC'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
22	vars(nv).fldName='DYC'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
23	vars(nv).fldName='DXG'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
24	vars(nv).fldName='DYG'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
25	vars(nv).fldName='hFacC'; vars(nv).longName='fractional thickness'; vars(nv).units='1'; nv=length(vars)+1;
26	vars(nv).fldName='hFacW'; vars(nv).longName='fractional thickness'; vars(nv).units='1'; nv=length(vars)+1;
27	vars(nv).fldName='hFacS'; vars(nv).longName='fractional thickness'; vars(nv).units='1'; nv=length(vars)+1;
28	vars(nv).fldName='Depth'; vars(nv).longName='sea floor depth'; vars(nv).units='m'; nv=length(vars)+1;
29	vars(nv).fldName='AngleCS'; vars(nv).longName='grid orientation (cosine)'; vars(nv).units='m'; nv=length(vars)+1;
30	vars(nv).fldName='AngleSN'; vars(nv).longName='grid orientation (sine)'; vars(nv).units='m'; nv=length(vars)+1;
31	vars(nv).fldName='RC'; vars(nv).longName='vertical coordinate'; vars(nv).units='m'; nv=length(vars)+1;
32	vars(nv).fldName='RF'; vars(nv).longName='vertical coordinate'; vars(nv).units='m'; nv=length(vars)+1;
33	vars(nv).fldName='DRC'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
34	vars(nv).fldName='DRF'; vars(nv).longName='grid spacing'; vars(nv).units='m'; nv=length(vars)+1;
35	%
36	structIn.defs=vars;
37	%
38	struct2nctiles('release1/','GRID',structIn,[90 90]);
39	end;
40
41	%===============
42
43	if choiceStruct==2;
44
45	listBudgs={'budgMo','budgHo','budgSo','budgMi','budgHi','budgSi'};
46
47	%directories and snapshot
48	dirIn='r4it11.c65i/';
49	dirMat='mat_budg3d/';
50	fldName='snapshot';
51	timName='0000000732'; suff='initial';
52	% timName='0000174564'; suff='final';
53	t0Name=num2str(3600*str2num(timName));
54
55	%load myparms
56	eval(['load ' dirIn dirMat 'diags_grid_parms.mat myparms;']);
57
58	%read THETA, SALT
59	diagName=['budg3d_snap_set1.' timName];
60	diagName=fullfile(dirIn,'diags',filesep,'BUDG',filesep,diagName);
61	THETA=rdmds2gcmfaces(diagName,'rec',1);
62	SALT=rdmds2gcmfaces(diagName,'rec',2);
63	ONE=1*(SALT>0);
64
65	%compute cell thickness
66	etanName=['budg2d_snap_set1.' timName];
67	etanName=fullfile(dirIn,'diags',filesep,'BUDG',filesep,etanName);
68	etan=rdmds2gcmfaces(etanName,'rec',1);
69	%compute time variable thickness
70	tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
71	tmp2=tmp1./mk3D(mygrid.Depth,tmp1);
72	tmp2=tmp2.*mk3D(etan,tmp1);
73	tmp3=mk3D(mygrid.RAC,tmp1);
74	vol=tmp3.*(tmp1+tmp2);
75
76	%load ice and snow volume/m2
77	SIheff=rdmds2gcmfaces(etanName,'rec',2);
78	SIsnow=rdmds2gcmfaces(etanName,'rec',3);
79
80	for ii=1:length(listBudgs);
81	budgName=listBudgs{ii};
82	[budgName ' -- ' fldName ' -- ' timName]
83
84	%directories and snapshot
85	dirOut=fullfile(dirIn,'nctiles_budg',filesep);
86	if ~isdir(dirOut); mkdir(dirOut); end;
87	dirOut=fullfile(dirOut,budgName,filesep);
88	if ~isdir(dirOut); mkdir(dirOut); end;
89
90	switch budgName;
91	case 'budgMo'; fld=ONE; fac=myparms.rhoconst;
92	case 'budgHo'; fld=THETA; fac=myparms.rcp;
93	case 'budgSo'; fld=SALT; fac=myparms.rhoconst;
94	case 'budgMi'; FACheff=myparms.rhoi; FACsnow=myparms.rhosn;
95	case 'budgHi'; FACheff=-myparms.flamimyparms.rhoi; FACsnow=-myparms.flamimyparms.rhosn;
96	case 'budgSi'; FACheff=myparms.SIsal0*myparms.rhoi; FACsnow=0;
97	end;
98
99	%compute extensive snapshot
100	structIn=[];
101	if budgName(end)=='o';
102	structIn.vars.snapshot=facfld.vol;
103	else;
104	structIn.vars.snapshot=mygrid.RAC.(FACheffSIheff+FACsnow*SIsnow);
105	end;
106
107	%general description
108	specs=[];
109	if strcmp(budgName(5),'H');
110	structIn.descr={'Heat budget in extensive form (in Watt, on C-Grid)'};
111	specs.units='J';
112	specs.name='heat content';
113	elseif strcmp(budgName(5),'M');
114	structIn.descr={'Mass budget in extensive form (in kg/s, on C-Grid)'};
115	specs.units='kg';
116	specs.name='mass content';
117	elseif strcmp(budgName(5),'S');
118	structIn.descr={'Salt budget in extensive form (in g/s, on C-Grid)'};
119	specs.units='g';
120	specs.name='salt content';
121	end;
122
123	%variables description
124	vars=[]; nv=length(vars)+1;
125	vars(nv).fldName=fldName; vars(nv).units=specs.units;
126	if strcmp(budgName(6),'o'); vars(nv).longName=['ocean ' specs.name ' snapshot'];
127	elseif strcmp(budgName(6),'i'); vars(nv).longName=['ice+snow ' specs.name ' snapshot'];
128	end;
129	vars(nv).longName=[vars(nv).longName ' at ' suff ' time=' t0Name 's'];
130	%
131	structIn.defs=vars;
132
133	%create file
134	tic; struct2nctiles(dirIn,fldName,structIn,[90 90]); toc;
135	eval(['!mv ' pwd filesep dirIn filesep 'tmp_nctiles' filesep fldName ' ' pwd filesep dirOut suff]);
136
137	end;%for ii=1:length(listBudgs);
138
139	end;
140
141	%===============
142
143	if choiceStruct==3;
144
145	%budgName='budgHo';
146
147	if budgName(end)=='o'; listFlds={'tend','trU','trV','trWtop'};
148	else; listFlds={'tend','trU','trV','trWtop','trWbot'};
149	end;
150
151	for ii=1:length(listFlds);
152	fldName=listFlds{ii};
153	%fldName='trWtop';
154	[budgName ' -- ' fldName]
155
156	%directories
157	dirIn='r4it11.c65i/';
158	dirMat='mat_budg3d/';
159	dirOut=fullfile(dirIn,'nctiles_budg',filesep);
160	if ~isdir(dirOut); mkdir(dirOut); end;
161	dirOut=fullfile(dirOut,budgName,filesep);
162	if ~isdir(dirOut); mkdir(dirOut); end;
163
164	%load variable
165	eval(['load ' dirIn dirMat 'diags_grid_parms.mat myparms;']);
166	dirMat=[dirIn dirMat 'diags_set_' budgName '/'];
167	fileMat=[budgName '_*.mat'];
168	tic; structIn.vars=diags_read_from_mat(dirMat,fileMat,fldName); toc;
169
170	%time vectors
171	[listTimes]=diags_list_times({[dirIn 'diags/BUDG/']},{'budg2d_hflux_set1'});
172	structIn.vars.t0=3600*listTimes(1:end-2);
173	structIn.vars.t1=3600*listTimes(2:end-1);
174
175	%rename trWtop as trW if adequate
176	if strcmp(fldName,'trWtop')&(budgName(end)=='o');
177	%note: since geothermal heating was added we likely need to add a nr+1 level to trW
178	% using trWbot(:,:,nr) or output trWtop and trWbot themselves (as done for seaice)
179	structIn.vars=setfield(structIn.vars,'trW',structIn.vars.trWtop);
180	structIn.vars=rmfield(structIn.vars,'trWtop');
181	structIn.vars.listDiags={'trW'};
182	fldName='trW';
183	end;
184
185	%switch back to upward convention
186	if strcmp(fldName(1:3),'trW');
187	%note: the following was needed in preparing the 03-Feb-2015 nctiles_budget (3d) files
188	% but is no longer needed since the trW sign convention in gcmfaces_diags was
189	% reversed to positive upward for consistency with MITgcm convention.
190	tmp1=getfield(structIn.vars,fldName);
191	structIn.vars=setfield(structIn.vars,fldName,-tmp1);
192	end;
193
194	%general description
195	tmp1=diags_read_from_mat(dirMat,fileMat,'specs',1);
196	specs=tmp1.specs;
197	if strcmp(specs.units,'W');
198	structIn.descr={'Heat budget in extensive form (in Watt, on C-Grid)'};
199	elseif strcmp(specs.units,'kg/s');
200	structIn.descr={'Mass budget in extensive form (in kg/s, on C-Grid)'};
201	elseif strcmp(specs.units,'g/s');
202	structIn.descr={'Salt budget in extensive form (in g/s, on C-Grid)'};
203	else;
204	error('unknown budget');
205	end;
206
207	%variables description
208	vars=[]; nv=length(vars)+1;
209	vars(nv).fldName=fldName; vars(nv).units=specs.units;
210	switch fldName;
211	case 'tend'; vars(nv).longName='tendency term';
212	case 'trU'; vars(nv).longName='horizontal transport (U)';
213	case 'trV'; vars(nv).longName='horizontal transport (V)';
214	case 'trW'; vars(nv).longName='upward vertical transport (W)';
215	case 'trWtop'; vars(nv).longName='upward vertical transport (W)';
216	case 'trWbot'; vars(nv).longName='upward vertical transport (W)';
217	end;
218	nv=length(vars)+1;
219	vars(nv).fldName='t0'; vars(nv).longName='initial time'; vars(nv).units='s'; nv=length(vars)+1;
220	vars(nv).fldName='t1'; vars(nv).longName='final time'; vars(nv).units='s'; nv=length(vars)+1;
221	%
222	structIn.defs=vars;
223
224	%create file
225	tic; struct2nctiles(dirIn,fldName,structIn,[90 90]); toc;
226	eval(['!mv ' pwd filesep dirIn filesep 'tmp_nctiles' filesep fldName ' ' pwd filesep dirOut]);
227
228	end;%for fldName=listFlds;
229
230	end;
231