clear all, close all

% domain-specific preamble
face=3; ix=326:365; jx=311:350; kx=1:50;   % domain specification
nme='beaufort';                            % domain name
obcs_src='cube78';                         % source of open boundaries
nt=179;                                    % number of obcs time steps

% directory names (may need to be created or modified)
eval(['cd /skylla/' nme]);
pout=['/skylla/' nme '/run_template/'];    % output path name
pin='/skylla/cube/run_template/';          % CS510 input files
grid='/skylla/cube/grid/cube66/';          % location of CS510 grid files
pnm=['/skylla/cube/' obcs_src '/'];        % open boundary path name

% miscellaneous input files
bathy_file='GEBCO_510x6x510_ver06_dig.bin';
diffkr_file= ...
    '/skylla/data/atmos/blend_forcing/cube78_forcing/DIFFKR_2_20_1_lat6070_cube78';

% derived quantities
nx=length(ix); ny=length(jx); nz=length(kx);
dim=[num2str(nx) 'x' num2str(ny)];

% CS510 grid input files
rc=-readbin([grid 'RC.data'],nz);          % depths to center of cell
rf=-readbin([grid 'RF.data'],nz+1);        % depths to cell faces
thk=diff(rf);                              % CS510 thicknesses
xc=read_cs_ifjk([grid 'XC.data'],ix,face,jx);
yc=read_cs_ifjk([grid 'YC.data'],ix,face,jx);

% location of domain-specific grid files: these need to be
% generated by running the model for at least one time step
% prior to balancing BCs, the last item in this script.
regional_grid=['/skylla/' nme '/grid/'];
genBC={'N','W','S'};                    % generate genBC boundary conditions
balBC='W';                              % balance balBC boundary condition


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create tile input
fn=[pin 'tile00' int2str(face) '.mitgrid'];
tile=readbin(fn,[511 511 16],1,'real*8');
writebin([pout 'LONC.bin'],tile(ix,jx, 1));
writebin([pout 'LATC.bin'],tile(ix,jx, 2));
writebin([pout 'DXF.bin' ],tile(ix,jx, 3));
writebin([pout 'DYF.bin' ],tile(ix,jx, 4));
writebin([pout 'RA.bin'  ],tile(ix,jx, 5));
writebin([pout 'LONG.bin'],tile(ix,jx, 6));
writebin([pout 'LATG.bin'],tile(ix,jx, 7));
writebin([pout 'DXV.bin' ],tile(ix,jx, 8));
writebin([pout 'DYU.bin' ],tile(ix,jx, 9));
writebin([pout 'RAZ.bin' ],tile(ix,jx,10));
writebin([pout 'DXC.bin' ],tile(ix,jx,11));
writebin([pout 'DYC.bin' ],tile(ix,jx,12));
writebin([pout 'RAW.bin' ],tile(ix,jx,13));
writebin([pout 'RAS.bin' ],tile(ix,jx,14));
writebin([pout 'DXG.bin' ],tile(ix,jx,15));
writebin([pout 'DYG.bin' ],tile(ix,jx,16));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make bathymetry file
topog=read_cs_ifjk([pin bathy_file],ix,face,jx);
if strcmp(nme,'arctic')
    topog(368:420,1:43)=0; topog(375:417,1:63)=0; topog(376:396,1:80)=0;
    topog(377:386,1:92)=0; topog(355:420,220:384)=0; topog(408:420,213:384)=0;
    topog(413:420,207:384)=0; topog(417:420,202:384)=0; topog(94:103,379:384)=0;
elseif strcmp(nme,'weddell')
    topog(37:41,1:7)=0; topog(48:54,1:12)=0;
end
writebin([pout 'BATHY_' obcs_src '_' dim '_' nme],topog);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make DIFFKR file
diffkr=read_cs_face(diffkr_file,face,kx);
fout=[pout 'DIFFKR_' obcs_src '_' dim 'x' int2str(length(kx)) '_' nme];
writebin(fout,diffkr(ix,jx,:));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make pickup files
% requires that pickup*.meta files be constructed mannually
fn=[pin 'pickup.0000000216.' obcs_src];
for k=1:303, mydisp(k)
 tmp=read_cs_face(fn,face,k,'real*8');
 writebin([pout 'pickup.0000000216.data'],tmp(ix,jx),1,'real*8',k-1);
end
fn=[pin 'pickup_seaice.0000000216.' obcs_src];
for k=1:22, mydisp(k)
 tmp=read_cs_face(fn,face,k,'real*8');
 writebin([pout 'pickup_seaice.0000000216.data'],tmp(ix,jx),1,'real*8',k-1);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make start-up files
for clm={'PHC','WGHC','WOA01','WOA05'}
    for fld={'SALT','THETA'}
        fin=[pin clm{1} '_' fld{1} '_JAN_510x6x510x50.bin'];
        fout=[pout clm{1} '_' fld{1} '_JAN_' dim  'x50_' nme];
        tmp=read_cs_face(fin,face,kx);
        writebin(fout,tmp(ix,jx,:));
    end
end
fn=[pin 'pickup_seaice.0000000216.' obcs_src];
id=[9 16 19 22]; fld={'HSNOW','HEFF','AREA','HSALT'}
for i=1:length(id)
 tmp=read_cs_face(fn,face,id(i),'real*8'); tmp(find(tmp<0))=0;
 if id(i)==19, tmp(find(tmp>1))=1; end
 writebin([pout fld{i} '_' dim '_' nme '.' obcs_src],tmp(ix,jx));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate seaice boundary conditions
for fld={'SIheff','SIarea','SIhsalt','SIhsnow','SIuice','SIvice'}
  switch fld{1}
    case 'SIheff',  id = 'h'   ;
    case 'SIarea',  id = 'a'   ;
    case 'SIhsalt', id = 'sl'  ;
    case 'SIhsnow', id = 'sn'  ;
    case 'SIuice',  id = 'uice';
    case 'SIvice',  id = 'vice';
  end
  fnm=dir([pnm fld{1} '/' fld{1} '.0*']);
  for i=1:length(fnm), mydisp(i)
    tmp=read_cs_face([pnm fld{1} '/' fnm(i).name],face);
    for g=1:length(genBC)
      switch genBC{g}
        case 'N', OB=tmp(ix,max(jx),:);
        case 'S'
          if strcmp(id,'vice')
            OB=tmp(ix,min(jx)+1,:);
          else
            OB=tmp(ix,min(jx),:);
          end
        case 'E', OB=tmp(max(ix),jx,:);
        case 'W'
          if strcmp(id,'uice')
            OB=tmp(min(ix)+1,jx,:);
          else
            OB=tmp(min(ix),jx,:);
          end
      end
      fout=[pout 'OB' genBC{g} id '_' nme '_' dim '.bin'];
      if i==1
        % add 4 identical daily records so that it is possible,
        % if desired, to start from beginning of 1992
        for r=1:4, writebin(fout,OB), end
      end
      writebin(fout,OB,1,'real*4',i+3)
    end
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate U/V/T/S lateral boundary conditions
for fld={'THETA','SALTanom','UVELMASS','VVELMASS'}
  fnm=dir([pnm fld{1} '/' fld{1} '.0*']);
  switch fld{1}
    case{'THETA','SALTanom'}, prec='real*8';
    otherwise prec='real*4';
  end
  for i=1:length(fnm), mydisp(i)
    tmp=read_cs_face([pnm fld{1} '/' fnm(i).name],face,kx,prec);
    if strcmp(fld{1},'SALTanom'), tmp=tmp+35; end
    for g=1:length(genBC)
      switch genBC{g}
        case 'N', OB=tmp(ix,max(jx),:);
        case 'S'
          if strcmp(fld{1},'VVELMASS')
            OB=tmp(ix,min(jx)+1,:);
          else
            OB=tmp(ix,min(jx),:);
          end
        case 'E', OB=tmp(max(ix),jx,:);
        case 'W'
          if strcmp(fld{1},'UVELMASS')
            OB=tmp(min(ix)+1,jx,:);
          else
            OB=tmp(min(ix),jx,:);
          end
      end
      fout=[pout 'OB' genBC{g} lower(fld{1}(1)) '_' nme '_' dim '.bin'];
      if i==1, writebin(fout,OB), end
      writebin(fout,OB,1,'real*4',i)
    end
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stabilize T/S
tmp=read_cs_face([grid 'hFacC.data'],face,1:50);
maskW=squeeze(tmp(ix(1),jx,:));
maskW(find(maskW))=1;  maskW(find(~maskW))=nan;
maskE=squeeze(tmp(max(ix),jx,:));
maskE(find(maskE))=1;  maskE(find(~maskE))=nan;
maskN=squeeze(tmp(ix,max(jx),:));
maskN(find(maskN))=1;  maskN(find(~maskN))=nan;
maskS=squeeze(tmp(ix,jx(1),:));
maskS(find(maskS))=1;  maskS(find(~maskS))=nan;
for t=1:nt, mydisp(t)
  for g=1:length(genBC)
    switch genBC{g}
      case 'N'
        T=readbin([pout 'OBNt_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1).*maskN;
        S=readbin([pout 'OBNs_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1).*maskN;
      case 'S'
        T=readbin([pout 'OBSt_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1).*maskS;
        S=readbin([pout 'OBSs_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1).*maskS;
      case 'E'
        T=readbin([pout 'OBEt_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1).*maskE;
        S=readbin([pout 'OBEs_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1).*maskE;
      case 'W'
        T=readbin([pout 'OBWt_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1).*maskW;
        S=readbin([pout 'OBWs_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1).*maskW;
    end
    R=rho(S,T,0);
    for j=1:ny
      idx=find(diff(R(j,:))<0);
      while ~isempty(idx)
        T(j,min(idx)+1)=T(j,min(idx));
        S(j,min(idx)+1)=S(j,min(idx));
        r=rho(S(j,:),T(j,:),0); idx=find(diff(r)<0);
      end
    end
    for k=1:nz
      if any(~isnan(T(:,k)))
        T(:,k)=xpolate(T(:,k)); S(:,k)=xpolate(S(:,k));
      else
        T(:,k)=T(:,k-1); S(:,k)=S(:,k-1);
      end
    end
    writebin([pout 'OB' genBC{g} 's_' nme '_' dim '.stable'],S,1,'real*4',t-1);
    writebin([pout 'OB' genBC{g} 't_' nme '_' dim '.stable'],T,1,'real*4',t-1);
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% balance BCs -- requires that model be integrated once so that
% grid and mask files be available.
DXG  =readbin([regional_grid 'DXG.data'   ],[nx ny]   );
DYG  =readbin([regional_grid 'DYG.data'   ],[nx ny]   );
HFACS=readbin([regional_grid 'hFacS.data'],[nx ny nz]);
HFACW=readbin([regional_grid 'hFacW.data'],[nx ny nz]);
OBWmask=squeeze(HFACW(2 ,: ,:)); OBWmask([1 ny],:)=0;
OBEmask=squeeze(HFACW(nx,: ,:)); OBEmask([1 ny],:)=0;
OBSmask=squeeze(HFACS(: ,2 ,:)); OBSmask([1 nx],:)=0;
OBNmask=squeeze(HFACS(: ,ny,:)); OBNmask([1 nx],:)=0;
for k=1:nz
  OBWmask(:,k)=thk(k)*OBWmask(:,k).*DYG(2 ,: )';
  OBEmask(:,k)=thk(k)*OBEmask(:,k).*DYG(nx,: )';
  OBSmask(:,k)=thk(k)*OBSmask(:,k).*DXG(: ,2 ) ;
  OBNmask(:,k)=thk(k)*OBNmask(:,k).*DXG(: ,ny) ;
end
OBW=zeros(nt,1); OBE=OBW; OBS=OBW; OBN=OBW;
for t=1:nt, mydisp(t)
    for g=1:length(genBC)
        switch genBC{g}
          case 'E'
            tmp=readbin([pout 'OBEu_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1);
            OBE(t)=sum(sum(tmp.*OBEmask));
          case 'W'
            tmp=readbin([pout 'OBWu_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1);
            OBW(t)=sum(sum(tmp.*OBWmask));
          case 'N'
            tmp=readbin([pout 'OBNv_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1);
            OBN(t)=sum(sum(tmp.*OBNmask));
          case 'S'
            tmp=readbin([pout 'OBSv_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1);
            OBS(t)=sum(sum(tmp.*OBSmask));
        end
    end
end
for t=1:nt, mydisp(t)
    switch balBC
      case 'W'
        tmp=readbin([pout 'OBWu_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1);
        tmp(find(OBWmask))=tmp(find(OBWmask))+ ...
            mean(OBE(t)-OBS(t)+OBN(t)-OBW(t))/sum(sum(OBWmask));
        writebin([pout 'OBWu_' nme '_' dim '.balance'],tmp,1,'real*4',t-1);
      case 'E'
        tmp=readbin([pout 'OBEu_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1);
        tmp(find(OBEmask))=tmp(find(OBEmask))+ ...
            mean(OBW(t)-OBN(t)+OBS(t)-OBE(t))/sum(sum(OBEmask));
        writebin([pout 'OBEu_' nme '_' dim '.balance'],tmp,1,'real*4',t-1);
      case 'N'
        tmp=readbin([pout 'OBNv_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1);
        tmp(find(OBNmask))=tmp(find(OBNmask))+ ...
            mean(-OBE(t)+OBW(t)+OBS(t)-OBN(t))/sum(sum(OBNmask));
        writebin([pout 'OBNv_' nme '_' dim '.balance'],tmp,1,'real*4',t-1);
      case 'S'
        tmp=readbin([pout 'OBSv_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1);
        tmp(find(OBSmask))=tmp(find(OBSmask))+ ...
            mean(OBE(t)-OBW(t)+OBN(t)-OBS(t))/sum(sum(OBSmask));
        writebin([pout 'OBSv_' nme '_' dim '.balance'],tmp,1,'real*4',t-1);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cd /skylla/beaufort/run_template
for t=1:nt, mydisp(t)
  tmp=readbin('OBNv_beaufort_40x40.bin',[nx nz],1,'real*4',t-1);
  OBN(t)=sum(sum(tmp.*OBNmask));
  tmp=readbin('OBSv_beaufort_40x40.bin',[nx nz],1,'real*4',t-1);
  OBS(t)=sum(sum(tmp.*OBSmask));
  tmp=readbin('OBWu_beaufort_40x40.bin',[ny nz],1,'real*4',t-1);
  OBW(t)=sum(sum(tmp.*OBWmask));
  tmp=readbin('OBWu_beaufort_40x40.balance',[ny nz],1,'real*4',t-1);
  OBW2(t)=sum(sum(tmp.*OBWmask));
end
t=(1:nt)';
clf, plot(t,cumsum(OBW2'-OBN+OBS),t,cumsum(OBW-OBN+OBS)), grid
     legend('cumsum(OBW-OBN2+OBS)','cumsum(OBW-OBN+OBS)')