Warning, /verification/so_box_biogeo/inp_global/mk_box_input.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit a1d0e455 on 2022-12-09 03:00:25 UTC
b94725a104 Jean*0001 % Generate input files for MITgcm/verification/so_box_biogeo/
0002 % ( Southern-Ocean box with Biochemistry, using OBCs )
0003 % by a) extracting bathy and forcing from ../run_glob/ where files have been
0004 % linked from verification/tutorial_global_oce_biogeo/input/
0005 % b) extracting initial conditions and OBCs from a 2-years run
0006 % (done here in ../run_glob/) using model parameters from ../inp_global.
0007
0008 %- kwr(1) = 2 : write Region mask to get matching Stat-Diags in Global Ocean run
0009 %- kwr(2) = 2 : write forcing files for S.Ocean box
0010 %- kwr(3) = 2 : write initial State for S.Ocean box
0011 %- kwr(4) = 2 : write OBCS files for S.Ocean box
0012 kwr=[1 2 2 2]; prec='real*4';
0013
0014 rDir='../run_glob/'; %- output dir of global ocean 2 yrs run
0015 iDir='../run_glob/'; %- input dir of global ocean 2 yrs run
0016
0017 G=load_grid(rDir);
0018 nx=G.dims(1); ny=G.dims(2); nr=G.dims(3);
0019 xax=G.xAxC; yax=G.yAxC;
0020 xaf=G.xAxU; yaf=G.yAxV;
0021
0022 it0=5184000;
0023 itMon=60; %- nb of iter in 1 month
0024
0025 nbx=42; nby=20;
0026 i1=85; i2=i1+nbx; %- box egdes indices
0027 j1=4; j2=j1+nby; %- box egdes indices
0028 xl=[xaf(i1) xaf(i1) xaf(i2) xaf(i2) xaf(i1)];
0029 yl=[yaf(j1) yaf(j2) yaf(j2) yaf(j1) yaf(j1)];
0030 ib1=i1; ib2=i2-1; jb1=j1; jb2=j2-1; %- box limit for grid-pts center
0031
0032 %-- OB location: tracer pts:
0033 iobW=ib1; iobE=ib2; jobN=jb2;
0034 %-- OB location: velocity, normal component:
0035 iuW=iobW+1 ; iuE=iobE ; jvN=jobN ;
0036 xob=[xaf(iuW) xaf(iuW) xaf(iuE) xaf(iuE)];
0037 yob=[yaf(j1) yaf(jvN) yaf(jvN) yaf(j1) ];
0038
0039 h=G.depth;
0040
0041 %- Create a Region mask (corresponding to S.Ocean Box interior)
0042 % to get matching Stat-Diags in Global Ocean run
0043 if kwr(1) > 0,
0044 % regional mask: 3 lat. band (=1,2,3) (y-boundary=-/+ 23) + S.Ocean Box (=0):
0045 maskDiag=ones(nx,ny);
0046 maskDiag(iobW+1:iobE-1,jb1:jobN-1)=0.;
0047 maskDiag(find(G.yC > -23.))=2.;
0048 maskDiag(find(G.yC > +23.))=3.;
0049 namfb='regMask_SOceanBox.bin';
0050 if kwr(1) > 1,
0051 vv2=maskDiag;
0052 fprintf(' write file: %-25s (%i %i)',namfb,size(vv2));
0053 fid=fopen(namfb,'w','b'); fwrite(fid,var,prec); fclose(fid);
0054 fprintf('\n');
0055 end
0056 end %- if kwr(1) > 0
0057
0058 figure(1);clf;
0059 var=-h; var(find(var==0))=NaN; ccB=[-54 1]*100;
0060 if kwr(1) > 0,
0061 var=maskDiag; var(find(h==0))=NaN; ccB=[-1 4];
0062 end
0063 imagesc(xax,yax,var');set(gca,'YDir','normal');
0064 caxis(ccB);
0065 change_colmap(-1);
0066 colorbar
0067 Lbx=line(xl,yl); set(Lbx,'Color',[0 0 0],'LineWidth',2);
0068 Lob=line(xob,yob); set(Lob,'Color',[1 0 0],'LineWidth',2);
0069 grid;
0070
0071 %---------
0072 %- extract box part from 2-D & 3-D fields
0073 % a) forcing (2-D bin file, sng or 12 Rec):
0074 namInpF={'bathy','lev_monthly_salt','lev_monthly_temp', ...
0075 'shi_empmr_year','shi_qnet','tren_taux','tren_tauy', ...
a1d0e455fd Hann*0076 'tren_speed','fice','sillev1','silicate_3D_12m'};
b94725a104 Jean*0077 if kwr(2) > 0
0078 fprintf('==== Processing Forcing Files ====\n');
0079
0080 for n=1:length(namInpF),
0081 namfg=[char(namInpF(n)),'.bin'];
0082 namfb=[char(namInpF(n)),'_box.bin'];
a1d0e455fd Hann*0083 if isempty(dir([iDir,namfg])),
0084 fprintf(' => skip file: %-21s (not found)\n',namfg);
0085 else
0086 fprintf(' process file: %-21s',namfg);
0087 nRec=12;
0088 if strcmp(char(namInpF(n)),'bathy'), nRec=1; end
0089 if strcmp(char(namInpF(n)),'silicate_3D_12m'), nRec=nr*nRec; end
0090 vv1=rdda([iDir,namfg],[nx ny nRec],1,prec,'b');
0091 vv2=vv1(ib1:ib2,jb1:jb2,:);
0092 %size(vv1),size(vv2)
0093 if kwr(2) > 1,
0094 fprintf(' --> %-25s (',namfb); fprintf(' %i',size(vv2)); fprintf(' )');
b94725a104 Jean*0095 fid=fopen(namfb,'w','b'); fwrite(fid,vv2,prec); fclose(fid);
a1d0e455fd Hann*0096 end
0097 fprintf('\n');
b94725a104 Jean*0098 end
0099 end
0100 end %- if kwr(2) > 0
0101
0102 % b) initial state
0103 % take snap-shot after 1 yr run:
0104 it1=it0+12*itMon;
0105 namState={'Eta','T','S','U','V'};
0106 if kwr(3) > 0
0107 fprintf('==== Processing Initial State ====\n');
0108
0109 for n=1:length(namState)+5,
0110 if n <= length(namState),
0111 namfg=[char(namState(n))];
0112 namfb=[char(namState(n)),'_ini.bin'];
0113 else
0114 itr= n-length(namState);
0115 namfg=sprintf('%s%2.2i','PTRACER',itr);
0116 namfb=sprintf('%s%2.2i%s','Trac',itr,'_ini.bin');
0117 end
0118 fprintf(' process file: %-21s',namfg);
0119 vv1=rdmds([rDir,namfg],it1);
0120 if length(size(vv1)) == 2,
0121 vv2=vv1(ib1:ib2,jb1:jb2);
0122 else
0123 vv2=vv1(ib1:ib2,jb1:jb2,:);
0124 end
0125 if kwr(3) > 1,
0126 fprintf(' --> %-25s',namfb);
0127 fid=fopen(namfb,'w','b'); fwrite(fid,vv2,prec); fclose(fid);
0128 end
0129 fprintf('\n');
0130 end
0131 end %- if kwr(3) > 0
0132
0133 %---------
0134 %- extract OB fields from 3-D output fields
0135 % take 1rst yr Dec aver as 1rst Rec ; then continue on to 2nd yr: Jan -> Nov
0136 iters=it0+([0:11]+12)*itMon;
0137
0138 if kwr(4) > 0,
0139 fprintf('==== Processing OB Fields ====\n');
0140 [vv4,itrs,M]=rdmds([rDir,'dynDiag'],iters);
0141 eval(M);
0142
0143 Nit=length(iters);
0144 rhFacW=G.hFacW; rhFacS=G.hFacS;
0145 rhFacW(find(rhFacW==0))=-1; rhFacW=1./rhFacW; rhFacW=max(rhFacW,0.);
0146 rhFacS(find(rhFacS==0))=-1; rhFacS=1./rhFacS; rhFacS=max(rhFacS,0.);
0147 rhFacW=reshape(rhFacW,[nx*ny*nr 1]);
0148 rhFacS=reshape(rhFacS,[nx*ny*nr 1]);
0149
0150 for n=2:length(namState)+5,
0151 if n > length(namState),
0152 itr= n-length(namState);
0153 namVs=sprintf('%s%2.2i','Trac',itr);
0154 namV=sprintf('%s%2.2i%s','TRAC',itr,' ');
0155 else
0156 namVs =[char(namState(n))];
0157 if strcmp(namVs,'U') || strcmp(namVs,'V'),
0158 namV=[namVs,'VELMASS'];
0159 elseif strcmp(namVs,'T'),
0160 namV=[namVs,'HETA '];
0161 elseif strcmp(namVs,'S'),
0162 namV=[namVs,'ALT '];
0163 %namV='WVELMASS';
0164 else
0165 fprintf(' unknow State Var: %s\n',namVs);
0166 end
0167 end
0168 jv=strmatch(namV,fldList);
0169 if isempty(jv),
0170 fprintf('not found in fldList: "%s" for Var: "%s"',namV,namVs);
0171 else
0172 fprintf(' process OB for: %-6s from %-8s (jv=%2i)',namVs,char(fldList(jv)),jv);
0173 vv1=squeeze(vv4(:,:,:,jv,:));
0174 if strcmp(namVs,'U'),
0175 vv1=reshape(vv1,[nx*ny*nr Nit]);
0176 vv1=vv1.*(rhFacW*ones(1,Nit));
0177 vv1=reshape(vv1,[nx ny nr Nit]);
0178 %-- extract OB values:
0179 vvW=vv1(iuW,jb1:jb2,:,:);
0180 vvE=vv1(iuE,jb1:jb2,:,:);
0181 vvN=vv1(ib1:ib2,jobN,:,:);
0182 elseif strcmp(namVs,'V'),
0183 vv1=reshape(vv1,[nx*ny*nr Nit]);
0184 vv1=vv1.*(rhFacS*ones(1,Nit));
0185 vv1=reshape(vv1,[nx ny nr Nit]);
0186 %-- extract OB values:
0187 vvW=vv1(iobW,jb1:jb2,:,:);
0188 vvE=vv1(iobE,jb1:jb2,:,:);
0189 vvN=vv1(ib1:ib2,jvN,:,:);
0190 else
0191 vv2=vv1(ib1:ib2,jb1:jb2,:);
0192 %-- extract Western OB values:
0193 vvW=vv1(iobW,jb1:jb2,:,:);
0194 %-- extract Eastern OB values:
0195 vvE=vv1(iobE,jb1:jb2,:,:);
0196 %-- extract Nortern OB values:
0197 vvN=vv1(ib1:ib2,jobN,:,:);
0198 end
0199 %--
0200 end
0201 if kwr(4) > 1,
0202 namfb=[namVs,'_obW.bin'];
0203 fprintf(' --> %s',namfb);
0204 fid=fopen(namfb,'w','b'); fwrite(fid,vvW,prec); fclose(fid);
0205 namfb=[namVs,'_obE.bin'];
0206 fprintf(' , %s',namfb);
0207 fid=fopen(namfb,'w','b'); fwrite(fid,vvE,prec); fclose(fid);
0208 namfb=[namVs,'_obN.bin'];
0209 fprintf(' , %s',namfb);
0210 fid=fopen(namfb,'w','b'); fwrite(fid,vvN,prec); fclose(fid);
0211 end
0212 fprintf('\n');
0213 end
0214 end %- if kwr(4) > 0
0215
