| 1 |
clear all, close all |
| 2 |
|
| 3 |
% domain-specific preamble |
| 4 |
face=3; ix=326:365; jx=311:350; kx=1:50; % domain specification |
| 5 |
nme='beaufort'; % domain name |
| 6 |
obcs_src='cube78'; % source of open boundaries |
| 7 |
nt=179; % number of obcs time steps |
| 8 |
|
| 9 |
% directory names (may need to be created or modified) |
| 10 |
eval(['cd /skylla/' nme]); |
| 11 |
pout=['/skylla/' nme '/run_template/']; % output path name |
| 12 |
pin='/skylla/cube/run_template/'; % CS510 input files |
| 13 |
grid='/skylla/cube/grid/cube66/'; % location of CS510 grid files |
| 14 |
pnm=['/skylla/cube/' obcs_src '/']; % open boundary path name |
| 15 |
|
| 16 |
% miscellaneous input files |
| 17 |
bathy_file='GEBCO_510x6x510_ver06_dig.bin'; |
| 18 |
diffkr_file= ... |
| 19 |
'/skylla/data/atmos/blend_forcing/cube78_forcing/DIFFKR_2_20_1_lat6070_cube78'; |
| 20 |
|
| 21 |
% derived quantities |
| 22 |
nx=length(ix); ny=length(jx); nz=length(kx); |
| 23 |
dim=[num2str(nx) 'x' num2str(ny)]; |
| 24 |
|
| 25 |
% CS510 grid input files |
| 26 |
rc=-readbin([grid 'RC.data'],nz); % depths to center of cell |
| 27 |
rf=-readbin([grid 'RF.data'],nz+1); % depths to cell faces |
| 28 |
thk=diff(rf); % CS510 thicknesses |
| 29 |
xc=read_cs_ifjk([grid 'XC.data'],ix,face,jx); |
| 30 |
yc=read_cs_ifjk([grid 'YC.data'],ix,face,jx); |
| 31 |
|
| 32 |
% location of domain-specific grid files: these need to be |
| 33 |
% generated by running the model for at least one time step |
| 34 |
% prior to balancing BCs, the last item in this script. |
| 35 |
regional_grid=['/skylla/' nme '/grid/']; |
| 36 |
genBC={'N','W','S'}; % generate genBC boundary conditions |
| 37 |
balBC='W'; % balance balBC boundary condition |
| 38 |
|
| 39 |
|
| 40 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 41 |
% create tile input |
| 42 |
fn=[pin 'tile00' int2str(face) '.mitgrid']; |
| 43 |
tile=readbin(fn,[511 511 16],1,'real*8'); |
| 44 |
writebin([pout 'LONC.bin'],tile(ix,jx, 1)); |
| 45 |
writebin([pout 'LATC.bin'],tile(ix,jx, 2)); |
| 46 |
writebin([pout 'DXF.bin' ],tile(ix,jx, 3)); |
| 47 |
writebin([pout 'DYF.bin' ],tile(ix,jx, 4)); |
| 48 |
writebin([pout 'RA.bin' ],tile(ix,jx, 5)); |
| 49 |
writebin([pout 'LONG.bin'],tile(ix,jx, 6)); |
| 50 |
writebin([pout 'LATG.bin'],tile(ix,jx, 7)); |
| 51 |
writebin([pout 'DXV.bin' ],tile(ix,jx, 8)); |
| 52 |
writebin([pout 'DYU.bin' ],tile(ix,jx, 9)); |
| 53 |
writebin([pout 'RAZ.bin' ],tile(ix,jx,10)); |
| 54 |
writebin([pout 'DXC.bin' ],tile(ix,jx,11)); |
| 55 |
writebin([pout 'DYC.bin' ],tile(ix,jx,12)); |
| 56 |
writebin([pout 'RAW.bin' ],tile(ix,jx,13)); |
| 57 |
writebin([pout 'RAS.bin' ],tile(ix,jx,14)); |
| 58 |
writebin([pout 'DXG.bin' ],tile(ix,jx,15)); |
| 59 |
writebin([pout 'DYG.bin' ],tile(ix,jx,16)); |
| 60 |
|
| 61 |
|
| 62 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 63 |
% make bathymetry file |
| 64 |
topog=read_cs_ifjk([pin bathy_file],ix,face,jx); |
| 65 |
if strcmp(nme,'arctic') |
| 66 |
topog(368:420,1:43)=0; topog(375:417,1:63)=0; topog(376:396,1:80)=0; |
| 67 |
topog(377:386,1:92)=0; topog(355:420,220:384)=0; topog(408:420,213:384)=0; |
| 68 |
topog(413:420,207:384)=0; topog(417:420,202:384)=0; topog(94:103,379:384)=0; |
| 69 |
elseif strcmp(nme,'weddell') |
| 70 |
topog(37:41,1:7)=0; topog(48:54,1:12)=0; |
| 71 |
end |
| 72 |
writebin([pout 'BATHY_' obcs_src '_' dim '_' nme],topog); |
| 73 |
|
| 74 |
|
| 75 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 76 |
% make DIFFKR file |
| 77 |
diffkr=read_cs_face(diffkr_file,face,kx); |
| 78 |
fout=[pout 'DIFFKR_' obcs_src '_' dim 'x' int2str(length(kx)) '_' nme]; |
| 79 |
writebin(fout,diffkr(ix,jx,:)); |
| 80 |
|
| 81 |
|
| 82 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 83 |
% make pickup files |
| 84 |
% requires that pickup*.meta files be constructed mannually |
| 85 |
fn=[pin 'pickup.0000000216.' obcs_src]; |
| 86 |
for k=1:303, mydisp(k) |
| 87 |
tmp=read_cs_face(fn,face,k,'real*8'); |
| 88 |
writebin([pout 'pickup.0000000216.data'],tmp(ix,jx),1,'real*8',k-1); |
| 89 |
end |
| 90 |
fn=[pin 'pickup_seaice.0000000216.' obcs_src]; |
| 91 |
for k=1:22, mydisp(k) |
| 92 |
tmp=read_cs_face(fn,face,k,'real*8'); |
| 93 |
writebin([pout 'pickup_seaice.0000000216.data'],tmp(ix,jx),1,'real*8',k-1); |
| 94 |
end |
| 95 |
|
| 96 |
|
| 97 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 98 |
% make start-up files |
| 99 |
for clm={'PHC','WGHC','WOA01','WOA05'} |
| 100 |
for fld={'SALT','THETA'} |
| 101 |
fin=[pin clm{1} '_' fld{1} '_JAN_510x6x510x50.bin']; |
| 102 |
fout=[pout clm{1} '_' fld{1} '_JAN_' dim 'x50_' nme]; |
| 103 |
tmp=read_cs_face(fin,face,kx); |
| 104 |
writebin(fout,tmp(ix,jx,:)); |
| 105 |
end |
| 106 |
end |
| 107 |
fn=[pin 'pickup_seaice.0000000216.' obcs_src]; |
| 108 |
id=[9 16 19 22]; fld={'HSNOW','HEFF','AREA','HSALT'} |
| 109 |
for i=1:length(id) |
| 110 |
tmp=read_cs_face(fn,face,id(i),'real*8'); tmp(find(tmp<0))=0; |
| 111 |
if id(i)==19, tmp(find(tmp>1))=1; end |
| 112 |
writebin([pout fld{i} '_' dim '_' nme '.' obcs_src],tmp(ix,jx)); |
| 113 |
end |
| 114 |
|
| 115 |
|
| 116 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 117 |
% generate seaice boundary conditions |
| 118 |
for fld={'SIheff','SIarea','SIhsalt','SIhsnow','SIuice','SIvice'} |
| 119 |
switch fld{1} |
| 120 |
case 'SIheff', id = 'h' ; |
| 121 |
case 'SIarea', id = 'a' ; |
| 122 |
case 'SIhsalt', id = 'sl' ; |
| 123 |
case 'SIhsnow', id = 'sn' ; |
| 124 |
case 'SIuice', id = 'uice'; |
| 125 |
case 'SIvice', id = 'vice'; |
| 126 |
end |
| 127 |
fnm=dir([pnm fld{1} '/' fld{1} '.0*']); |
| 128 |
for i=1:length(fnm), mydisp(i) |
| 129 |
tmp=read_cs_face([pnm fld{1} '/' fnm(i).name],face); |
| 130 |
for g=1:length(genBC) |
| 131 |
switch genBC{g} |
| 132 |
case 'N', OB=tmp(ix,max(jx),:); |
| 133 |
case 'S' |
| 134 |
if strcmp(id,'vice') |
| 135 |
OB=tmp(ix,min(jx)+1,:); |
| 136 |
else |
| 137 |
OB=tmp(ix,min(jx),:); |
| 138 |
end |
| 139 |
case 'E', OB=tmp(max(ix),jx,:); |
| 140 |
case 'W' |
| 141 |
if strcmp(id,'uice') |
| 142 |
OB=tmp(min(ix)+1,jx,:); |
| 143 |
else |
| 144 |
OB=tmp(min(ix),jx,:); |
| 145 |
end |
| 146 |
end |
| 147 |
fout=[pout 'OB' genBC{g} id '_' nme '_' dim '.bin']; |
| 148 |
if i==1 |
| 149 |
% add 4 identical daily records so that it is possible, |
| 150 |
% if desired, to start from beginning of 1992 |
| 151 |
for r=1:4, writebin(fout,OB), end |
| 152 |
end |
| 153 |
writebin(fout,OB,1,'real*4',i+3) |
| 154 |
end |
| 155 |
end |
| 156 |
end |
| 157 |
|
| 158 |
|
| 159 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 160 |
% generate U/V/T/S lateral boundary conditions |
| 161 |
for fld={'THETA','SALTanom','UVELMASS','VVELMASS'} |
| 162 |
fnm=dir([pnm fld{1} '/' fld{1} '.0*']); |
| 163 |
switch fld{1} |
| 164 |
case{'THETA','SALTanom'}, prec='real*8'; |
| 165 |
otherwise prec='real*4'; |
| 166 |
end |
| 167 |
for i=1:length(fnm), mydisp(i) |
| 168 |
tmp=read_cs_face([pnm fld{1} '/' fnm(i).name],face,kx,prec); |
| 169 |
if strcmp(fld{1},'SALTanom'), tmp=tmp+35; end |
| 170 |
for g=1:length(genBC) |
| 171 |
switch genBC{g} |
| 172 |
case 'N', OB=tmp(ix,max(jx),:); |
| 173 |
case 'S' |
| 174 |
if strcmp(fld{1},'VVELMASS') |
| 175 |
OB=tmp(ix,min(jx)+1,:); |
| 176 |
else |
| 177 |
OB=tmp(ix,min(jx),:); |
| 178 |
end |
| 179 |
case 'E', OB=tmp(max(ix),jx,:); |
| 180 |
case 'W' |
| 181 |
if strcmp(fld{1},'UVELMASS') |
| 182 |
OB=tmp(min(ix)+1,jx,:); |
| 183 |
else |
| 184 |
OB=tmp(min(ix),jx,:); |
| 185 |
end |
| 186 |
end |
| 187 |
fout=[pout 'OB' genBC{g} lower(fld{1}(1)) '_' nme '_' dim '.bin']; |
| 188 |
if i==1, writebin(fout,OB), end |
| 189 |
writebin(fout,OB,1,'real*4',i) |
| 190 |
end |
| 191 |
end |
| 192 |
end |
| 193 |
|
| 194 |
|
| 195 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 196 |
% stabilize T/S |
| 197 |
tmp=read_cs_face([grid 'hFacC.data'],face,1:50); |
| 198 |
maskW=squeeze(tmp(ix(1),jx,:)); |
| 199 |
maskW(find(maskW))=1; maskW(find(~maskW))=nan; |
| 200 |
maskE=squeeze(tmp(max(ix),jx,:)); |
| 201 |
maskE(find(maskE))=1; maskE(find(~maskE))=nan; |
| 202 |
maskN=squeeze(tmp(ix,max(jx),:)); |
| 203 |
maskN(find(maskN))=1; maskN(find(~maskN))=nan; |
| 204 |
maskS=squeeze(tmp(ix,jx(1),:)); |
| 205 |
maskS(find(maskS))=1; maskS(find(~maskS))=nan; |
| 206 |
for t=1:nt, mydisp(t) |
| 207 |
for g=1:length(genBC) |
| 208 |
switch genBC{g} |
| 209 |
case 'N' |
| 210 |
T=readbin([pout 'OBNt_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1).*maskN; |
| 211 |
S=readbin([pout 'OBNs_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1).*maskN; |
| 212 |
case 'S' |
| 213 |
T=readbin([pout 'OBSt_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1).*maskS; |
| 214 |
S=readbin([pout 'OBSs_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1).*maskS; |
| 215 |
case 'E' |
| 216 |
T=readbin([pout 'OBEt_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1).*maskE; |
| 217 |
S=readbin([pout 'OBEs_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1).*maskE; |
| 218 |
case 'W' |
| 219 |
T=readbin([pout 'OBWt_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1).*maskW; |
| 220 |
S=readbin([pout 'OBWs_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1).*maskW; |
| 221 |
end |
| 222 |
R=rho(S,T,0); |
| 223 |
for j=1:ny |
| 224 |
idx=find(diff(R(j,:))<0); |
| 225 |
while ~isempty(idx) |
| 226 |
T(j,min(idx)+1)=T(j,min(idx)); |
| 227 |
S(j,min(idx)+1)=S(j,min(idx)); |
| 228 |
r=rho(S(j,:),T(j,:),0); idx=find(diff(r)<0); |
| 229 |
end |
| 230 |
end |
| 231 |
for k=1:nz |
| 232 |
if any(~isnan(T(:,k))) |
| 233 |
T(:,k)=xpolate(T(:,k)); S(:,k)=xpolate(S(:,k)); |
| 234 |
else |
| 235 |
T(:,k)=T(:,k-1); S(:,k)=S(:,k-1); |
| 236 |
end |
| 237 |
end |
| 238 |
writebin([pout 'OB' genBC{g} 's_' nme '_' dim '.stable'],S,1,'real*4',t-1); |
| 239 |
writebin([pout 'OB' genBC{g} 't_' nme '_' dim '.stable'],T,1,'real*4',t-1); |
| 240 |
end |
| 241 |
end |
| 242 |
|
| 243 |
|
| 244 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 245 |
% balance BCs -- requires that model be integrated once so that |
| 246 |
% grid and mask files be available. |
| 247 |
DXG =readbin([regional_grid 'DXG.data' ],[nx ny] ); |
| 248 |
DYG =readbin([regional_grid 'DYG.data' ],[nx ny] ); |
| 249 |
HFACS=readbin([regional_grid 'hFacS.data'],[nx ny nz]); |
| 250 |
HFACW=readbin([regional_grid 'hFacW.data'],[nx ny nz]); |
| 251 |
OBWmask=squeeze(HFACW(2 ,: ,:)); OBWmask([1 ny],:)=0; |
| 252 |
OBEmask=squeeze(HFACW(nx,: ,:)); OBEmask([1 ny],:)=0; |
| 253 |
OBSmask=squeeze(HFACS(: ,2 ,:)); OBSmask([1 nx],:)=0; |
| 254 |
OBNmask=squeeze(HFACS(: ,ny,:)); OBNmask([1 nx],:)=0; |
| 255 |
for k=1:nz |
| 256 |
OBWmask(:,k)=thk(k)*OBWmask(:,k).*DYG(2 ,: )'; |
| 257 |
OBEmask(:,k)=thk(k)*OBEmask(:,k).*DYG(nx,: )'; |
| 258 |
OBSmask(:,k)=thk(k)*OBSmask(:,k).*DXG(: ,2 ) ; |
| 259 |
OBNmask(:,k)=thk(k)*OBNmask(:,k).*DXG(: ,ny) ; |
| 260 |
end |
| 261 |
OBW=zeros(nt,1); OBE=OBW; OBS=OBW; OBN=OBW; |
| 262 |
for t=1:nt, mydisp(t) |
| 263 |
for g=1:length(genBC) |
| 264 |
switch genBC{g} |
| 265 |
case 'E' |
| 266 |
tmp=readbin([pout 'OBEu_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1); |
| 267 |
OBE(t)=sum(sum(tmp.*OBEmask)); |
| 268 |
case 'W' |
| 269 |
tmp=readbin([pout 'OBWu_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1); |
| 270 |
OBW(t)=sum(sum(tmp.*OBWmask)); |
| 271 |
case 'N' |
| 272 |
tmp=readbin([pout 'OBNv_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1); |
| 273 |
OBN(t)=sum(sum(tmp.*OBNmask)); |
| 274 |
case 'S' |
| 275 |
tmp=readbin([pout 'OBSv_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1); |
| 276 |
OBS(t)=sum(sum(tmp.*OBSmask)); |
| 277 |
end |
| 278 |
end |
| 279 |
end |
| 280 |
for t=1:nt, mydisp(t) |
| 281 |
switch balBC |
| 282 |
case 'W' |
| 283 |
tmp=readbin([pout 'OBWu_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1); |
| 284 |
tmp(find(OBWmask))=tmp(find(OBWmask))+ ... |
| 285 |
mean(OBE(t)-OBS(t)+OBN(t)-OBW(t))/sum(sum(OBWmask)); |
| 286 |
writebin([pout 'OBWu_' nme '_' dim '.balance'],tmp,1,'real*4',t-1); |
| 287 |
case 'E' |
| 288 |
tmp=readbin([pout 'OBEu_' nme '_' dim '.bin'],[ny nz],1,'real*4',t-1); |
| 289 |
tmp(find(OBEmask))=tmp(find(OBEmask))+ ... |
| 290 |
mean(OBW(t)-OBN(t)+OBS(t)-OBE(t))/sum(sum(OBEmask)); |
| 291 |
writebin([pout 'OBEu_' nme '_' dim '.balance'],tmp,1,'real*4',t-1); |
| 292 |
case 'N' |
| 293 |
tmp=readbin([pout 'OBNv_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1); |
| 294 |
tmp(find(OBNmask))=tmp(find(OBNmask))+ ... |
| 295 |
mean(-OBE(t)+OBW(t)+OBS(t)-OBN(t))/sum(sum(OBNmask)); |
| 296 |
writebin([pout 'OBNv_' nme '_' dim '.balance'],tmp,1,'real*4',t-1); |
| 297 |
case 'S' |
| 298 |
tmp=readbin([pout 'OBSv_' nme '_' dim '.bin'],[nx nz],1,'real*4',t-1); |
| 299 |
tmp(find(OBSmask))=tmp(find(OBSmask))+ ... |
| 300 |
mean(OBE(t)-OBW(t)+OBN(t)-OBS(t))/sum(sum(OBSmask)); |
| 301 |
writebin([pout 'OBSv_' nme '_' dim '.balance'],tmp,1,'real*4',t-1); |
| 302 |
end |
| 303 |
end |
| 304 |
|
| 305 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 306 |
|
| 307 |
cd /skylla/beaufort/run_template |
| 308 |
for t=1:nt, mydisp(t) |
| 309 |
tmp=readbin('OBNv_beaufort_40x40.bin',[nx nz],1,'real*4',t-1); |
| 310 |
OBN(t)=sum(sum(tmp.*OBNmask)); |
| 311 |
tmp=readbin('OBSv_beaufort_40x40.bin',[nx nz],1,'real*4',t-1); |
| 312 |
OBS(t)=sum(sum(tmp.*OBSmask)); |
| 313 |
tmp=readbin('OBWu_beaufort_40x40.bin',[ny nz],1,'real*4',t-1); |
| 314 |
OBW(t)=sum(sum(tmp.*OBWmask)); |
| 315 |
tmp=readbin('OBWu_beaufort_40x40.balance',[ny nz],1,'real*4',t-1); |
| 316 |
OBW2(t)=sum(sum(tmp.*OBWmask)); |
| 317 |
end |
| 318 |
t=(1:nt)'; |
| 319 |
clf, plot(t,cumsum(OBW2'-OBN+OBS),t,cumsum(OBW-OBN+OBS)), grid |
| 320 |
legend('cumsum(OBW-OBN2+OBS)','cumsum(OBW-OBN+OBS)') |
Parent Directory
| 
Revision Log
| 
Revision Graph