% vgrid.m creates the model vertical grid and stores the grid information
% various arrays. Nominally set are:
%    level thicknesses (dz)
% from which are generated:
%    tracer  grid coordinates (xc,yc,zc)
%    w-point grid coordinates (xc,yc,zf)
%
%      load VGRID.mat                   (loads everything above)
%  
% Uses data files: n/a
% Creates data files: VGRID.mat
% Creates arrays: dz,zc,zf
%
% Created  04/20/01 by adcroft@mit.edu

% Vertical level thicknesses
% 30 levels, -6565
%dz=[10 10 12.5 15 17.5 20 25 30 35 40 50 60 70 80 100 120 140 170 200 240 280 330 370 420 480 540 600 660 700 740];
% 20 levels, -6000 (exp2)
%dz=[50 50 55 60 65 70 80 95 120 155 200 260 320 400 480 570 655 725 775 815];
% 15 levels, -6000 (global_ocean)
%dz=[50 70 100 140 190 240 290 340 390 440 490 540 590 640 690];
% 23 levels -5700 m global ocean ("old" ecco configuration)
dz=[10.,10.,15.,20.,20.,25.,35.,50.,75.,100.,150.,200.,275.,350.,415., ...
    450.,500.,500.,500.,500.,500.,500.,500.];
% ecco-godae/ecco
% $$$ dz=[10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
% $$$     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04, 19.82, 24.85, ...
% $$$     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
% $$$     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
% $$$     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
% $$$     341.50,364.50,387.50,410.50,433.50,456.50];

nzc=prod(size(dz));

zf=-cumsum([0 dz]);
zc=(zf(1:end-1)+zf(2:end))/2;
disp(sprintf('nzc = %i',nzc))
disp(['dz =' sprintf(' %4.1f',dz)])

% Store data in VGRID.mat
save VGRID.mat nzc dz zc zf