[MITgcm-devel] experiment to check conservation in seaice

[Jean-Michel Campin]

Hi Martin,

I've checked the set-up, looks OK to me.
(relatively simple on the ocean side since there is no
relaxation on SSST & SSS).
I've changed check_conserve.m (attached), to check over the period
here it's 1 year) instead of per unit of time (/s), 
since there was already some conversion to /year, it's simpler.

I am getting a good conservation for the ocean (vol,salt,heat),
but did not look at seaice & atmos interface.

Cheers,
Jean-Michel

On Tue, Nov 02, 2010 at 03:23:13PM +0100, Martin Losch wrote:
> Hi Jean-Michel,
> 
> I finally set up a simple experiment based on the verification/global_ocean.cs32x15 with seaice and seaice thermodynamics (no winton). I failed miserably in trying to check the salt and heat conservation even for the ocean. The only thing that i managed to get right is volume conservation in the ocean, which is almost trivial of course. The seaice model has a drift of about 4mm/y as expected from earlier experiment.
> 
> I put the configuration on faulks:~mlosch/cs32conserve.tgz (just the code and input directories).
> 
> In the input directory there's a matlab script that tries to make use of the diagnostics output in netcdf (check_conserve.m). Please teach me how to close the heat and salt budgets.
> 
> Martin
> 
> 
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
-------------- next part --------------
bdir = '../run';
bdir = 'res_mnc';

% load the data first
if 0
spval = 123456.7;
si = rdmnc(fullfile(bdir,'SIceStDiag.*'));
ss = rdmnc(fullfile(bdir,'surfStDiag.*'));
is = rdmnc(fullfile(bdir,'instStDiag.*'));
fldsn=fieldnames(ss);
for k=1:length(fldsn)
  if ~strcmp(fldsn{k},'attributes');
    tmp = getfield(ss,fldsn{k});
    tmp(tmp==spval)=0;
    ss = setfield(ss,fldsn{k},tmp);
  end
end
fldsn=fieldnames(si);
for k=1:length(fldsn)
  if ~strcmp(fldsn{k},'attributes');
    tmp = getfield(si,fldsn{k});
    tmp(tmp==spval)=0;
    si = setfield(si,fldsn{k},tmp);
  end
end
fldsn=fieldnames(is);
for k=1:length(fldsn)
  if ~strcmp(fldsn{k},'attributes');
    tmp = getfield(is,fldsn{k});
    tmp(tmp==spval)=0;
    is = setfield(is,fldsn{k},tmp);
  end
end
end

% some constants
deltat= 360*86400;
rhosw = 1035; % sea water (=rhoConst)
rhofw = 1000;
rhoi  = 910;
rhosn = 330;
ru2m  = 1035;
rcp = rhosw*3994; % rho*(heat capacity)

% volume/mass
glbArea=is.ETAN_vol(1);
delEta = squeeze(   is.ETAN_ave(1,:,2)-   is.ETAN_ave(1,:,1));
delEta = squeeze(   is.ETAN_ave(1,:,2)-   is.ETAN_ave(1,:,1));
delHi  = squeeze( is.SIheff_ave(1,:,2)- is.SIheff_ave(1,:,1));
delHs  = squeeze(is.SIhsnow_ave(1,:,2)-is.SIhsnow_ave(1,:,1));

units=1; ustr = 'kg/m^2';
units=1000/1000; ustr = 'mm';
disp('global volume budget: ocean')
disp(sprintf('  rhosw*dEta = %u, oceFWflx*dt = %u', ...
	     rhosw*delEta(1)*units,ss.oceFWflx_ave(1)*deltat*units))
disp(sprintf('  residual = %u %s', ...
	     (rhosw*delEta(1)-ss.oceFWflx_ave(1)*deltat)*units,ustr))
disp('global volume budget: sea ice')
disp(sprintf('  rhoi*dHi+rhos*dHs = %u, net fwFlx *dt = %u', ...
	     (rhoi*delHi(1)+rhosn*delHs(1))*units, ...
	     (si.SIatmFW_ave(1)+si.SIempmr_ave(1))*deltat*units ) )
disp(sprintf('  residual = %u %s', ...
	     ((rhoi*delHi(1)+rhosn*delHs(1)) ...
	      -(si.SIatmFW_ave(1)+si.SIempmr_ave(1))*deltat)*units, ustr ...
	     )     )

disp('global volume budget: ocean + sea ice')
disp(sprintf('  rhosw*dEta+rhoi*dHi+rhos*dHs = %u, atm fwFlx *dt= %u',...
	     (rhosw*delEta(1)+rhoi*delHi(1)+rhosn*delHs(1))*units, ...
	     si.SIatmFW_ave(1)*deltat*units))
disp(sprintf('  residual = %u %s', ...
	     (si.SIatmFW_ave(1)*deltat ...
	     - (rhosw*delEta(1)+rhoi*delHi(1)+rhosn*delHs(1)))*units, ...
	     ustr ...
	     ))

% salt

units = 1e-9; ustr = 'ktons';
r3d = is.SALT_vol(1)*units;
r2d = ss.oceFWflx_vol(1)*units;
dsdt = squeeze(   is.SALT_ave(1,:,2)-   is.SALT_ave(1,:,1))/deltat;
disp('global salt budget: ocean')
disp(sprintf('  ru2m*dS/dt = %u, surForcS = %u, SFLUX = %u',...
	     ru2m*dsdt(1)*r3d, ss.surForcS_ave(1)*r2d, ss.SFLUX_ave(1)*r2d ) )
disp(sprintf('  residual = %u %s', ...
	     (ru2m*dsdt(1)-ss.surForcS_ave(1))*r2d,ustr) )

units = 1; ustr = 'g/m^2';
r3d=rhosw/glbArea;
dSalt = squeeze( is.SALT_ave(1,:,2).*is.SALT_vol(1,:,2) ...
               - is.SALT_ave(1,:,1).*is.SALT_vol(1,:,1) );
disp('-New- salt budget: ocean')
disp(sprintf('  dSalt = %u, oceSflux = %u, SFLUX = %u',...
	     dSalt(1)*r3d*units, ss.oceSflux_ave(1)*deltat*units, ...
                                 ss.SFLUX_ave(1)*deltat*units ) )
disp(sprintf('  residual = %u %s', ...
	     (dSalt(1)*r3d-ss.SFLUX_ave(1)*deltat)*units,ustr) )

% heat
units = 1e-12; ustr = 'TW';
dtdt = squeeze(   is.THETA_ave(1,:,2)-   is.THETA_ave(1,:,1))/deltat;
r3d = is.THETA_vol(1)*units;
r2d = ss.oceQnet_vol(1)*units;
disp('global heat budget: ocean')
disp(sprintf('  rcp*dT/dt = %u, oceQnet= %u, TFLUX = %u',...
	     rcp*dtdt(1)*r3d, ss.oceQnet_ave(1)*r2d, ...
	     ss.TFLUX_ave(1)*r2d ) )
disp(sprintf('  residual = %u %s', ...
	     (rcp*dtdt(1)*r3d - ss.oceQnet_ave(1)*r2d), ...
	     ustr ) )

units = 1; ustr = 'J/m^2';
r3d=rcp/glbArea;
dTheta = squeeze( is.THETA_ave(1,:,2).*is.THETA_vol(1,:,2) ...
                - is.THETA_ave(1,:,1).*is.THETA_vol(1,:,1) );
disp('-New- heat budget: ocean')
disp(sprintf('  rcp*dT = %u, oceQnet= %u, TFLUX = %u',...
	     dTheta(1)*r3d*units, ss.oceQnet_ave(1)*deltat*units, ...
	     ss.TFLUX_ave(1)*deltat*units ) )
disp(sprintf('  residual = %u %s', ...
	     (dTheta(1)*r3d - ss.oceQnet_ave(1)*deltat)*units, ...
	     ustr ) )

return