so_box_biogeo/inp_global/mk_box_input.m

% Generate input files for MITgcm/verification/so_box_biogeo/
%  ( Southern-Ocean box with Biochemistry, using OBCs )
% by a) extracting bathy and forcing from ../run_glob/ where files have been
%       linked from verification/tutorial_global_oce_biogeo/input/
%    b) extracting initial conditions and OBCs from a 2-years run
%       (done here in ../run_glob/) using model parameters from ../inp_global.

% $Header:  $
% $Name:  $

%- kwr(1) = 2 : write Region mask to get matching Stat-Diags in Global Ocean run
%- kwr(2) = 2 : write forcing files for S.Ocean box
%- kwr(3) = 2 : write initial State for S.Ocean box
%- kwr(4) = 2 : write OBCS files for S.Ocean box
kwr=[1 2 2 2]; prec='real*4';

rDir='../run_glob/'; %- output dir of global ocean 2 yrs run
iDir='../run_glob/'; %- input  dir of global ocean 2 yrs run

G=load_grid(rDir);
nx=G.dims(1); ny=G.dims(2); nr=G.dims(3);
xax=G.xAxC; yax=G.yAxC;
xaf=G.xAxU; yaf=G.yAxV;

it0=5184000;
itMon=60; %- nb of iter in 1 month

nbx=42; nby=20;
i1=85; i2=i1+nbx; %- box egdes indices
j1=4;  j2=j1+nby; %- box egdes indices
xl=[xaf(i1) xaf(i1) xaf(i2) xaf(i2) xaf(i1)];
yl=[yaf(j1) yaf(j2) yaf(j2) yaf(j1) yaf(j1)];
ib1=i1; ib2=i2-1; jb1=j1; jb2=j2-1; %- box limit for grid-pts center

%-- OB location: tracer pts:
iobW=ib1; iobE=ib2; jobN=jb2;
%-- OB location: velocity, normal component:
iuW=iobW+1 ; iuE=iobE ; jvN=jobN ;
xob=[xaf(iuW) xaf(iuW) xaf(iuE) xaf(iuE)];
yob=[yaf(j1)  yaf(jvN) yaf(jvN) yaf(j1) ];

h=G.depth;

%- Create a Region mask (corresponding to S.Ocean Box interior)
%    to get matching Stat-Diags in Global Ocean run
if kwr(1) > 0,
% regional mask: 3 lat. band (=1,2,3) (y-boundary=-/+ 23) + S.Ocean Box (=0):
  maskDiag=ones(nx,ny);
  maskDiag(iobW+1:iobE-1,jb1:jobN-1)=0.;
  maskDiag(find(G.yC > -23.))=2.;
  maskDiag(find(G.yC > +23.))=3.;
  namfb='regMask_SOceanBox.bin';
  if kwr(1) > 1,
   vv2=maskDiag;
   fprintf(' write file: %-25s (%i %i)',namfb,size(vv2));
   fid=fopen(namfb,'w','b'); fwrite(fid,var,prec); fclose(fid);
   fprintf('\n');
  end
end %- if kwr(1) > 0

figure(1);clf;
var=-h; var(find(var==0))=NaN; ccB=[-54 1]*100;
if kwr(1) > 0,
  var=maskDiag; var(find(h==0))=NaN; ccB=[-1 4];
end
imagesc(xax,yax,var');set(gca,'YDir','normal');
caxis(ccB);
change_colmap(-1);
colorbar
Lbx=line(xl,yl); set(Lbx,'Color',[0 0 0],'LineWidth',2);
Lob=line(xob,yob); set(Lob,'Color',[1 0 0],'LineWidth',2);
grid;

%---------
%- extract box part from 2-D & 3-D fields
% a) forcing (2-D bin file, sng or 12 Rec):
namInpF={'bathy','lev_monthly_salt','lev_monthly_temp', ...
         'shi_empmr_year','shi_qnet','tren_taux','tren_tauy', ...
         'tren_speed','fice','sillev1'};
if kwr(2) > 0
 fprintf('==== Processing Forcing Files ====\n');

 for n=1:length(namInpF),
  namfg=[char(namInpF(n)),'.bin'];
  namfb=[char(namInpF(n)),'_box.bin'];
  fprintf(' process file: %-21s',namfg);
  nRec=12; if n ==1, nRec=1; end
  vv1=rdda([iDir,namfg],[nx ny nRec],1,prec,'b');
  vv2=vv1(ib1:ib2,jb1:jb2,:);
 %size(vv1),size(vv2)
  if kwr(2) > 1,
    if nRec == 1,
      fprintf(' --> %-25s (%i %i)',namfb,size(vv2));
    else
      fprintf(' --> %-25s (%i %i %i)',namfb,size(vv2));
    end
    fid=fopen(namfb,'w','b'); fwrite(fid,vv2,prec); fclose(fid);
  end
  fprintf('\n');
 end
end %- if kwr(2) > 0

% b) initial state
% take snap-shot after 1 yr run:
it1=it0+12*itMon;
namState={'Eta','T','S','U','V'};
if kwr(3) > 0
 fprintf('==== Processing Initial State ====\n');

 for n=1:length(namState)+5,
  if n <= length(namState),
   namfg=[char(namState(n))];
   namfb=[char(namState(n)),'_ini.bin'];
  else
   itr= n-length(namState);
   namfg=sprintf('%s%2.2i','PTRACER',itr);
   namfb=sprintf('%s%2.2i%s','Trac',itr,'_ini.bin');
  end
  fprintf(' process file: %-21s',namfg);
  vv1=rdmds([rDir,namfg],it1);
  if length(size(vv1)) == 2,
   vv2=vv1(ib1:ib2,jb1:jb2);
  else
   vv2=vv1(ib1:ib2,jb1:jb2,:);
  end
  if kwr(3) > 1,
   fprintf(' --> %-25s',namfb);
   fid=fopen(namfb,'w','b'); fwrite(fid,vv2,prec); fclose(fid);
  end
  fprintf('\n');
 end
end %- if kwr(3) > 0

%---------
%- extract OB fields from 3-D output fields
%  take 1rst yr Dec aver as 1rst Rec ; then continue on to 2nd yr: Jan -> Nov
iters=it0+([0:11]+12)*itMon;

if kwr(4) > 0,
 fprintf('==== Processing OB Fields ====\n');
 [vv4,itrs,M]=rdmds([rDir,'dynDiag'],iters);
 eval(M);

 Nit=length(iters);
 rhFacW=G.hFacW; rhFacS=G.hFacS;
 rhFacW(find(rhFacW==0))=-1; rhFacW=1./rhFacW; rhFacW=max(rhFacW,0.);
 rhFacS(find(rhFacS==0))=-1; rhFacS=1./rhFacS; rhFacS=max(rhFacS,0.);
 rhFacW=reshape(rhFacW,[nx*ny*nr 1]);
 rhFacS=reshape(rhFacS,[nx*ny*nr 1]);

 for n=2:length(namState)+5,
  if n > length(namState),
   itr= n-length(namState);
   namVs=sprintf('%s%2.2i','Trac',itr);
   namV=sprintf('%s%2.2i%s','TRAC',itr,'  ');
  else
   namVs =[char(namState(n))];
   if strcmp(namVs,'U') || strcmp(namVs,'V'),
     namV=[namVs,'VELMASS'];
   elseif strcmp(namVs,'T'),
     namV=[namVs,'HETA   '];
   elseif strcmp(namVs,'S'),
     namV=[namVs,'ALT    '];
    %namV='WVELMASS';
   else
     fprintf(' unknow State Var: %s\n',namVs);
   end
  end
  jv=strmatch(namV,fldList);
  if isempty(jv),
   fprintf('not found in fldList: "%s" for Var: "%s"',namV,namVs);
  else
   fprintf(' process OB for: %-6s from %-8s (jv=%2i)',namVs,char(fldList(jv)),jv);
   vv1=squeeze(vv4(:,:,:,jv,:));
   if strcmp(namVs,'U'),
     vv1=reshape(vv1,[nx*ny*nr Nit]);
     vv1=vv1.*(rhFacW*ones(1,Nit));
     vv1=reshape(vv1,[nx ny nr Nit]);
%-- extract OB values:
     vvW=vv1(iuW,jb1:jb2,:,:);
     vvE=vv1(iuE,jb1:jb2,:,:);
     vvN=vv1(ib1:ib2,jobN,:,:);
   elseif strcmp(namVs,'V'),
     vv1=reshape(vv1,[nx*ny*nr Nit]);
     vv1=vv1.*(rhFacS*ones(1,Nit));
     vv1=reshape(vv1,[nx ny nr Nit]);
%-- extract OB values:
     vvW=vv1(iobW,jb1:jb2,:,:);
     vvE=vv1(iobE,jb1:jb2,:,:);
     vvN=vv1(ib1:ib2,jvN,:,:);
   else
   vv2=vv1(ib1:ib2,jb1:jb2,:);
%-- extract Western OB values:
     vvW=vv1(iobW,jb1:jb2,:,:);
%-- extract Eastern OB values:
     vvE=vv1(iobE,jb1:jb2,:,:);
%-- extract Nortern OB values:
     vvN=vv1(ib1:ib2,jobN,:,:);
   end
%--
  end
  if kwr(4) > 1,
   namfb=[namVs,'_obW.bin'];
   fprintf(' --> %s',namfb);
   fid=fopen(namfb,'w','b'); fwrite(fid,vvW,prec); fclose(fid);
   namfb=[namVs,'_obE.bin'];
   fprintf(' , %s',namfb);
   fid=fopen(namfb,'w','b'); fwrite(fid,vvE,prec); fclose(fid);
   namfb=[namVs,'_obN.bin'];
   fprintf(' , %s',namfb);
   fid=fopen(namfb,'w','b'); fwrite(fid,vvN,prec); fclose(fid);
  end
  fprintf('\n');
 end
end %- if kwr(4) > 0

1	jmc	1.1	% Generate input files for MITgcm/verification/so_box_biogeo/
2			% ( Southern-Ocean box with Biochemistry, using OBCs )
3			% by a) extracting bathy and forcing from ../run_glob/ where files have been
4			% linked from verification/tutorial_global_oce_biogeo/input/
5			% b) extracting initial conditions and OBCs from a 2-years run
6			% (done here in ../run_glob/) using model parameters from ../inp_global.
7
8			% $Header: $
9			% $Name: $
10
11			%- kwr(1) = 2 : write Region mask to get matching Stat-Diags in Global Ocean run
12			%- kwr(2) = 2 : write forcing files for S.Ocean box
13			%- kwr(3) = 2 : write initial State for S.Ocean box
14			%- kwr(4) = 2 : write OBCS files for S.Ocean box
15			kwr=[1 2 2 2]; prec='real*4';
16
17			rDir='../run_glob/'; %- output dir of global ocean 2 yrs run
18			iDir='../run_glob/'; %- input dir of global ocean 2 yrs run
19
20			G=load_grid(rDir);
21			nx=G.dims(1); ny=G.dims(2); nr=G.dims(3);
22			xax=G.xAxC; yax=G.yAxC;
23			xaf=G.xAxU; yaf=G.yAxV;
24
25			it0=5184000;
26			itMon=60; %- nb of iter in 1 month
27
28			nbx=42; nby=20;
29			i1=85; i2=i1+nbx; %- box egdes indices
30			j1=4; j2=j1+nby; %- box egdes indices
31			xl=[xaf(i1) xaf(i1) xaf(i2) xaf(i2) xaf(i1)];
32			yl=[yaf(j1) yaf(j2) yaf(j2) yaf(j1) yaf(j1)];
33			ib1=i1; ib2=i2-1; jb1=j1; jb2=j2-1; %- box limit for grid-pts center
34
35			%-- OB location: tracer pts:
36			iobW=ib1; iobE=ib2; jobN=jb2;
37			%-- OB location: velocity, normal component:
38			iuW=iobW+1 ; iuE=iobE ; jvN=jobN ;
39			xob=[xaf(iuW) xaf(iuW) xaf(iuE) xaf(iuE)];
40			yob=[yaf(j1) yaf(jvN) yaf(jvN) yaf(j1) ];
41
42			h=G.depth;
43
44			%- Create a Region mask (corresponding to S.Ocean Box interior)
45			% to get matching Stat-Diags in Global Ocean run
46			if kwr(1) > 0,
47			% regional mask: 3 lat. band (=1,2,3) (y-boundary=-/+ 23) + S.Ocean Box (=0):
48			maskDiag=ones(nx,ny);
49			maskDiag(iobW+1:iobE-1,jb1:jobN-1)=0.;
50			maskDiag(find(G.yC > -23.))=2.;
51			maskDiag(find(G.yC > +23.))=3.;
52			namfb='regMask_SOceanBox.bin';
53			if kwr(1) > 1,
54			vv2=maskDiag;
55			fprintf(' write file: %-25s (%i %i)',namfb,size(vv2));
56			fid=fopen(namfb,'w','b'); fwrite(fid,var,prec); fclose(fid);
57			fprintf('\n');
58			end
59			end %- if kwr(1) > 0
60
61			figure(1);clf;
62			var=-h; var(find(var==0))=NaN; ccB=[-54 1]*100;
63			if kwr(1) > 0,
64			var=maskDiag; var(find(h==0))=NaN; ccB=[-1 4];
65			end
66			imagesc(xax,yax,var');set(gca,'YDir','normal');
67			caxis(ccB);
68			change_colmap(-1);
69			colorbar
70			Lbx=line(xl,yl); set(Lbx,'Color',[0 0 0],'LineWidth',2);
71			Lob=line(xob,yob); set(Lob,'Color',[1 0 0],'LineWidth',2);
72			grid;
73
74			%---------
75			%- extract box part from 2-D & 3-D fields
76			% a) forcing (2-D bin file, sng or 12 Rec):
77			namInpF={'bathy','lev_monthly_salt','lev_monthly_temp', ...
78			'shi_empmr_year','shi_qnet','tren_taux','tren_tauy', ...
79			'tren_speed','fice','sillev1'};
80			if kwr(2) > 0
81			fprintf('==== Processing Forcing Files ====\n');
82
83			for n=1:length(namInpF),
84			namfg=[char(namInpF(n)),'.bin'];
85			namfb=[char(namInpF(n)),'_box.bin'];
86			fprintf(' process file: %-21s',namfg);
87			nRec=12; if n ==1, nRec=1; end
88			vv1=rdda([iDir,namfg],[nx ny nRec],1,prec,'b');
89			vv2=vv1(ib1:ib2,jb1:jb2,:);
90			%size(vv1),size(vv2)
91			if kwr(2) > 1,
92			if nRec == 1,
93			fprintf(' --> %-25s (%i %i)',namfb,size(vv2));
94			else
95			fprintf(' --> %-25s (%i %i %i)',namfb,size(vv2));
96			end
97			fid=fopen(namfb,'w','b'); fwrite(fid,vv2,prec); fclose(fid);
98			end
99			fprintf('\n');
100			end
101			end %- if kwr(2) > 0
102
103			% b) initial state
104			% take snap-shot after 1 yr run:
105			it1=it0+12*itMon;
106			namState={'Eta','T','S','U','V'};
107			if kwr(3) > 0
108			fprintf('==== Processing Initial State ====\n');
109
110			for n=1:length(namState)+5,
111			if n <= length(namState),
112			namfg=[char(namState(n))];
113			namfb=[char(namState(n)),'_ini.bin'];
114			else
115			itr= n-length(namState);
116			namfg=sprintf('%s%2.2i','PTRACER',itr);
117			namfb=sprintf('%s%2.2i%s','Trac',itr,'_ini.bin');
118			end
119			fprintf(' process file: %-21s',namfg);
120			vv1=rdmds([rDir,namfg],it1);
121			if length(size(vv1)) == 2,
122			vv2=vv1(ib1:ib2,jb1:jb2);
123			else
124			vv2=vv1(ib1:ib2,jb1:jb2,:);
125			end
126			if kwr(3) > 1,
127			fprintf(' --> %-25s',namfb);
128			fid=fopen(namfb,'w','b'); fwrite(fid,vv2,prec); fclose(fid);
129			end
130			fprintf('\n');
131			end
132			end %- if kwr(3) > 0
133
134			%---------
135			%- extract OB fields from 3-D output fields
136			% take 1rst yr Dec aver as 1rst Rec ; then continue on to 2nd yr: Jan -> Nov
137			iters=it0+([0:11]+12)*itMon;
138
139			if kwr(4) > 0,
140			fprintf('==== Processing OB Fields ====\n');
141			[vv4,itrs,M]=rdmds([rDir,'dynDiag'],iters);
142			eval(M);
143
144			Nit=length(iters);
145			rhFacW=G.hFacW; rhFacS=G.hFacS;
146			rhFacW(find(rhFacW==0))=-1; rhFacW=1./rhFacW; rhFacW=max(rhFacW,0.);
147			rhFacS(find(rhFacS==0))=-1; rhFacS=1./rhFacS; rhFacS=max(rhFacS,0.);
148			rhFacW=reshape(rhFacW,[nxnynr 1]);
149			rhFacS=reshape(rhFacS,[nxnynr 1]);
150
151			for n=2:length(namState)+5,
152			if n > length(namState),
153			itr= n-length(namState);
154			namVs=sprintf('%s%2.2i','Trac',itr);
155			namV=sprintf('%s%2.2i%s','TRAC',itr,' ');
156			else
157			namVs =[char(namState(n))];
158			if strcmp(namVs,'U') \|\| strcmp(namVs,'V'),
159			namV=[namVs,'VELMASS'];
160			elseif strcmp(namVs,'T'),
161			namV=[namVs,'HETA '];
162			elseif strcmp(namVs,'S'),
163			namV=[namVs,'ALT '];
164			%namV='WVELMASS';
165			else
166			fprintf(' unknow State Var: %s\n',namVs);
167			end
168			end
169			jv=strmatch(namV,fldList);
170			if isempty(jv),
171			fprintf('not found in fldList: "%s" for Var: "%s"',namV,namVs);
172			else
173			fprintf(' process OB for: %-6s from %-8s (jv=%2i)',namVs,char(fldList(jv)),jv);
174			vv1=squeeze(vv4(:,:,:,jv,:));
175			if strcmp(namVs,'U'),
176			vv1=reshape(vv1,[nxnynr Nit]);
177			vv1=vv1.(rhFacWones(1,Nit));
178			vv1=reshape(vv1,[nx ny nr Nit]);
179			%-- extract OB values:
180			vvW=vv1(iuW,jb1:jb2,:,:);
181			vvE=vv1(iuE,jb1:jb2,:,:);
182			vvN=vv1(ib1:ib2,jobN,:,:);
183			elseif strcmp(namVs,'V'),
184			vv1=reshape(vv1,[nxnynr Nit]);
185			vv1=vv1.(rhFacSones(1,Nit));
186			vv1=reshape(vv1,[nx ny nr Nit]);
187			%-- extract OB values:
188			vvW=vv1(iobW,jb1:jb2,:,:);
189			vvE=vv1(iobE,jb1:jb2,:,:);
190			vvN=vv1(ib1:ib2,jvN,:,:);
191			else
192			vv2=vv1(ib1:ib2,jb1:jb2,:);
193			%-- extract Western OB values:
194			vvW=vv1(iobW,jb1:jb2,:,:);
195			%-- extract Eastern OB values:
196			vvE=vv1(iobE,jb1:jb2,:,:);
197			%-- extract Nortern OB values:
198			vvN=vv1(ib1:ib2,jobN,:,:);
199			end
200			%--
201			end
202			if kwr(4) > 1,
203			namfb=[namVs,'_obW.bin'];
204			fprintf(' --> %s',namfb);
205			fid=fopen(namfb,'w','b'); fwrite(fid,vvW,prec); fclose(fid);
206			namfb=[namVs,'_obE.bin'];
207			fprintf(' , %s',namfb);
208			fid=fopen(namfb,'w','b'); fwrite(fid,vvE,prec); fclose(fid);
209			namfb=[namVs,'_obN.bin'];
210			fprintf(' , %s',namfb);
211			fid=fopen(namfb,'w','b'); fwrite(fid,vvN,prec); fclose(fid);
212			end
213			fprintf('\n');
214			end
215			end %- if kwr(4) > 0
216