| 1 |
gforget |
1.1 |
|
| 2 |
|
|
%'budgMo' ocean mass budget (volume*rhoconst) |
| 3 |
|
|
%'budgHo' ocean heat budget |
| 4 |
|
|
%'budgSo' ocean salt budget |
| 5 |
|
|
%'budgMi' seaice+snow mass budget |
| 6 |
|
|
%'budgHi' seaice+snow heat budget |
| 7 |
|
|
%'budgSi' seaice+snow salt budget |
| 8 |
|
|
|
| 9 |
|
|
budgName='budgHo' |
| 10 |
|
|
nt=238 |
| 11 |
|
|
nk=10; |
| 12 |
|
|
doRes=0; |
| 13 |
|
|
|
| 14 |
|
|
dirIn=[dir0 '/release1/nctiles_budg/' budgName '/']; |
| 15 |
|
|
|
| 16 |
|
|
%constant parameters: |
| 17 |
|
|
parms.rhoconst =1029; %sea water density |
| 18 |
|
|
parms.rcp =3994*parms.rhoconst; % sea water rho X heat capacity |
| 19 |
|
|
parms.rhoi = 910; %sea ice density |
| 20 |
|
|
parms.rhosn = 330; %snow density |
| 21 |
|
|
parms.flami = 3.34e05; % latent heat of fusion of ice/snow (J/kg) |
| 22 |
|
|
parms.flamb = 2.50e06; % latent heat of evaporation (J/kg) |
| 23 |
|
|
|
| 24 |
|
|
%reference grid cell volumes: |
| 25 |
|
|
tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC; |
| 26 |
|
|
tmp3=mk3D(mygrid.RAC,tmp1); |
| 27 |
|
|
vol0=tmp3.*(tmp1); |
| 28 |
|
|
|
| 29 |
|
|
switch budgName; |
| 30 |
|
|
case 'budgMo'; fac=parms.rhoconst; |
| 31 |
|
|
case 'budgHo'; fac=parms.rcp; |
| 32 |
|
|
case 'budgSo'; fac=parms.rhoconst; |
| 33 |
|
|
case 'budgMi'; FACheff=myparms.rhoi; FACsnow=myparms.rhosn; |
| 34 |
|
|
case 'budgHi'; FACheff=-myparms.flami*myparms.rhoi; FACsnow=-myparms.flami*myparms.rhosn; |
| 35 |
|
|
case 'budgSi'; FACheff=myparms.SIsal0*myparms.rhoi; FACsnow=0; |
| 36 |
|
|
end; |
| 37 |
|
|
|
| 38 |
|
|
budgIn=[]; |
| 39 |
|
|
|
| 40 |
|
|
budgIn.ini=read_nctiles([dirIn 'initial/snapshot'],'snapshot'); |
| 41 |
|
|
budgIn.fin=read_nctiles([dirIn 'final/snapshot'],'snapshot'); |
| 42 |
|
|
eval(['ncload ' dirIn 'tend/tend.0001.nc t0']); |
| 43 |
|
|
eval(['ncload ' dirIn 'tend/tend.0001.nc t1']); |
| 44 |
|
|
budgIn.dt=t1-t0; |
| 45 |
|
|
|
| 46 |
|
|
increments=repmat(NaN*mygrid.RAC,[1 1 nt]); |
| 47 |
|
|
|
| 48 |
|
|
switch budgName; |
| 49 |
|
|
case 'budgMo'; |
| 50 |
|
|
ini=nansum(budgIn.ini-fac*vol0,3)./mygrid.RAC/fac; |
| 51 |
|
|
fin=nansum(budgIn.fin-fac*vol0,3)./mygrid.RAC/fac; |
| 52 |
|
|
case 'budgHo'; |
| 53 |
|
|
ini=nansum(budgIn.ini(:,:,1:nk),3)./nansum(vol0(:,:,1:nk),3)/fac; |
| 54 |
|
|
fin=nansum(budgIn.fin(:,:,1:nk),3)./nansum(vol0(:,:,1:nk),3)/fac; |
| 55 |
|
|
otherwise; error('not implemented yet'); |
| 56 |
|
|
end; |
| 57 |
|
|
|
| 58 |
|
|
for tt=1:nt; |
| 59 |
|
|
disp(tt) |
| 60 |
|
|
%load the various fields: |
| 61 |
|
|
budgIn.tend=read_nctiles([dirIn 'tend/tend'],'tend',tt); |
| 62 |
|
|
budgIn.trU=read_nctiles([dirIn 'trU/trU'],'trU',tt); |
| 63 |
|
|
budgIn.trV=read_nctiles([dirIn 'trV/trV'],'trV',tt); |
| 64 |
|
|
if dirIn(end-1)=='o'; |
| 65 |
|
|
budgIn.trW=read_nctiles([dirIn 'trW/trW'],'trW',tt); |
| 66 |
|
|
budgIn.trWtop=budgIn.trW; |
| 67 |
|
|
budgIn.trWbot=budgIn.trW(:,:,2:50); |
| 68 |
|
|
budgIn.trWbot(:,:,50)=0; |
| 69 |
|
|
else; |
| 70 |
|
|
budgIn.trWtop=read_nctiles([dirIn 'trWtop/trWtop'],'trWtop',tt); |
| 71 |
|
|
budgIn.trWbot=read_nctiles([dirIn 'trWbot/trWbot'],'trWbot',tt); |
| 72 |
|
|
end; |
| 73 |
|
|
%compute budget residuals: |
| 74 |
|
|
if doRes; |
| 75 |
|
|
nr=size(budgIn.tend{1},3); |
| 76 |
|
|
for kk=1:nr; prec(kk,:)=check_budg(budgIn,kk); end; |
| 77 |
|
|
if tt==1; |
| 78 |
|
|
store_prec=prec; |
| 79 |
|
|
else; |
| 80 |
|
|
store_prec(:,:,tt)=prec; |
| 81 |
|
|
end; |
| 82 |
|
|
end; |
| 83 |
|
|
%reconstruct time series from budget |
| 84 |
|
|
switch budgName; |
| 85 |
|
|
case 'budgMo'; |
| 86 |
|
|
increments(:,:,tt)=budgIn.dt(tt)*nansum(budgIn.tend,3)./mygrid.RAC/fac; |
| 87 |
|
|
case 'budgHo'; |
| 88 |
|
|
increments(:,:,tt)=budgIn.dt(tt)*nansum(budgIn.tend(:,:,1:nk),3)... |
| 89 |
|
|
./nansum(vol0(:,:,1:nk),3)/fac; |
| 90 |
|
|
otherwise; error('not implemented yet'); |
| 91 |
|
|
end; |
| 92 |
|
|
|
| 93 |
|
|
end;%for tt=[1 119 238]; |
| 94 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph