matlab_class/gcmfaces_devel/idma_plot_budgets.m


%'budgMo' ocean mass budget (volume*rhoconst)
%'budgHo' ocean heat budget
%'budgSo' ocean salt budget
%'budgMi' seaice+snow mass budget
%'budgHi' seaice+snow heat budget
%'budgSi' seaice+snow salt budget

budgName='budgHo'
nt=238
nk=10;
doRes=0;

dirIn=[dir0 '/release1/nctiles_budg/' budgName '/'];

%constant parameters:
parms.rhoconst =1029; %sea water density
parms.rcp      =3994*parms.rhoconst; % sea water rho X heat capacity
parms.rhoi     = 910; %sea ice density
parms.rhosn    = 330; %snow density
parms.flami    = 3.34e05; % latent heat of fusion of ice/snow (J/kg)
parms.flamb    = 2.50e06; % latent heat of evaporation (J/kg)

%reference grid cell volumes:
tmp1=mk3D(mygrid.DRF,mygrid.hFacC).*mygrid.hFacC;
tmp3=mk3D(mygrid.RAC,tmp1);
vol0=tmp3.*(tmp1);

switch budgName;
case 'budgMo'; fac=parms.rhoconst;
case 'budgHo'; fac=parms.rcp;
case 'budgSo'; fac=parms.rhoconst;
case 'budgMi'; FACheff=myparms.rhoi; FACsnow=myparms.rhosn;
case 'budgHi'; FACheff=-myparms.flami*myparms.rhoi; FACsnow=-myparms.flami*myparms.rhosn;
case 'budgSi'; FACheff=myparms.SIsal0*myparms.rhoi; FACsnow=0;
end;

budgIn=[];

budgIn.ini=read_nctiles([dirIn 'initial/snapshot'],'snapshot');
budgIn.fin=read_nctiles([dirIn 'final/snapshot'],'snapshot');
eval(['ncload ' dirIn 'tend/tend.0001.nc t0']);
eval(['ncload ' dirIn 'tend/tend.0001.nc t1']);
budgIn.dt=t1-t0;

increments=repmat(NaN*mygrid.RAC,[1 1 nt]);

switch budgName;
case 'budgMo'; 
  ini=nansum(budgIn.ini-fac*vol0,3)./mygrid.RAC/fac;
  fin=nansum(budgIn.fin-fac*vol0,3)./mygrid.RAC/fac;
case 'budgHo';
  ini=nansum(budgIn.ini(:,:,1:nk),3)./nansum(vol0(:,:,1:nk),3)/fac;
  fin=nansum(budgIn.fin(:,:,1:nk),3)./nansum(vol0(:,:,1:nk),3)/fac;
otherwise; error('not implemented yet');
end;

for tt=1:nt;
  disp(tt)
  %load the various fields:
  budgIn.tend=read_nctiles([dirIn 'tend/tend'],'tend',tt);
  budgIn.trU=read_nctiles([dirIn 'trU/trU'],'trU',tt);
  budgIn.trV=read_nctiles([dirIn 'trV/trV'],'trV',tt);
  if dirIn(end-1)=='o';
    budgIn.trW=read_nctiles([dirIn 'trW/trW'],'trW',tt);
    budgIn.trWtop=budgIn.trW;
    budgIn.trWbot=budgIn.trW(:,:,2:50);
    budgIn.trWbot(:,:,50)=0;
  else;
    budgIn.trWtop=read_nctiles([dirIn 'trWtop/trWtop'],'trWtop',tt);
    budgIn.trWbot=read_nctiles([dirIn 'trWbot/trWbot'],'trWbot',tt);
  end;
  %compute budget residuals:
  if doRes;
    nr=size(budgIn.tend{1},3);
    for kk=1:nr; prec(kk,:)=check_budg(budgIn,kk); end;
    if tt==1;
      store_prec=prec;
    else;
      store_prec(:,:,tt)=prec;
    end;
  end;
  %reconstruct time series from budget
  switch budgName;
  case 'budgMo'; 
    increments(:,:,tt)=budgIn.dt(tt)*nansum(budgIn.tend,3)./mygrid.RAC/fac;
  case 'budgHo';
    increments(:,:,tt)=budgIn.dt(tt)*nansum(budgIn.tend(:,:,1:nk),3)...
                       ./nansum(vol0(:,:,1:nk),3)/fac;
  otherwise; error('not implemented yet');
  end;

end;%for tt=[1 119 238];

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.flamimyparms.rhoi; FACsnow=-myparms.flamimyparms.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