Isothermal layers volume budget with OCCA


This is a set of matlab and fortran routines to compute a thermally defined layer volume budget over a subdomain of the global ocean from the OCCA state estimate.

This package was developed by Guillaume Maze and Gael Forget at MIT (Cambridge, USA) for the ECCO group and the CLIMODE project.

You can browse the package here:

(See below for downloading)


To perform this analysis you'll need:

These are normally defined on a 360x160x50 grid points, daily from 2003/11/1 to 2006/11/3 (1099 days). You can download OCCA here: or email Gael Forget.

svn checkout Tlayer_budget


There are 2 steps in the computation. Step 1 extracts a subdomain of the Atlas, Step 2 computes the budget with a 3D interpolation of all fields.

Step 1

Variables that need to be tuned are:

You may want to tune the subfunction wherearewe to match your machine. You may also want to use a special mask for your domain of analysis. You can do so by looking at the subfunction extract. Finally, if you have a fancy OCCA dataset, you may have to tune how the routine name the input files. Check at variables prefi, suffi and the subfunction open_files. Last, the function extract_subdomain.m takes as an argument a number which correspond to a particular set of fields to extract:

Step 2

Then you have to select parameters for klist and tlist variables. They select what are the field(s) you want to compute the budget of and the time lines to use.

The variable klist can contain one of the following:

and tlist can contain:

jpi = diff(subdomain(1:2))+1;

jpj = diff(subdomain(3:4))+1;

jpk = diff(subdomain(5:6))+1;

Interpolation parameters are nb_in and nbIn2Out. nb_in is the factor by which the program will multiply the number of grid points. Starting from the 1x1 degree Atlas, nb_in=12 will interpolate the field to a 1/12 deg grid. The other parameter nbIn2Out sets the factor by which the interpolated fields are scaled for the output fields. nbIn2Out=3 will produce records on a 12/3=4, ie 1/4 degree grid. IMPORTANT: both nb_in and the result of nb_in/nbIn2Out must be even numbers (multiple of 2).


This is a bunch of routines to look at the results in the subfolder Tlayer_budget/extra. If you want to do it on your own, you should know what the output files are.

For each field, you have two output files: the one with 'timeseries' in the name corresponds to the time serie of the domain integrated field, the other the 3D map of the field (time serie cumulated).

For example in Matlab, you can try:

% Your parameters here:

iter0 = 1;

iter1 = 1099;

nt = iter1 - iter0 + 1;

THETAlow = 16;

THETAhigh = 19;

nb_in = 12;

nbIn2Out = 3;

nb_out = nb_in/nbIn2Out;

dirtofiles = '~/my_outputs_are_here/';

% If you used the default set up of file names, this should be ok:

prefi = sprintf('%s/InterpALLonce_%4.4d_%4.4d_%2.0f.m%2.0f._timeseries_',dirtofiles,iter0,iter1,THETAlow,THETAhig);

sg = -1/3; % This is because fields are cumulated over 3 years and on the other side of the temperature tendency equation

t = datenum(2003,11,1,0,0,0); t = t(1):t(1)+nt-1; % Time axis

% Read outputs:

fid = fopen(strcat(prefi,'LRvol.3yV2adv.bin'),'r','ieee-be');

Vvol = fread(fid,[1 nt],'float32'); % Contains the layer volume time serie

fid = fopen(strcat(prefi,'LRvol.3yV2adv.bin'),'r','ieee-be');

M = fread(fid,[jpi*jpj*nb_out*nb_out jpk*nb_out],'float32');fclose(fid);

M = reshape(M,[jpi*nb_out jpj*nb_out jpk*nb_out]);

Mvol = permute(M*sg,[3 2 1]); clear M

% Mvol(x,y,z) contains the ''concentration'' of the layer

fid = fopen(strcat(prefi,'LRairsea3D.3yV2adv.bin'),'r','ieee-be');

Vq = fread(fid,[1 nt],'float32'); % Contains the air-sea heat flux formation rate time serie

fid = fopen(strcat(prefi,'LRairsea3D.3yV2adv.bin'),'r','ieee-be');

M = fread(fid,[jpi*jpj*nb_out*nb_out jpk*nb_out],'float32');fclose(fid);

M = reshape(M,[jpi*nb_out jpj*nb_out jpk*nb_out]);

Mq = permute(M*sg,[3 2 1]); clear M

% Mq contains the map of formation rate due to air-sea heat fluxes.

% Note that the total integral of Mq must match the last value of Vq.

% For all fields, the unit is Sv.y = 10^6*86400*365 ~ 3.15x10^13 m^3