Back to home page

MITgcm

 
 

    


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