Documentation of get_w_from_divu

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

  Load u and v data;

Cross-Reference Information

Listing of script get_w_from_divu

clear
cd /home/disk/hayes2/dvimont/csiro/data

tim = 201:300;  %  First check to see if the routine works
lev = 1:9;
lims = [-0.1 360 -90 90];

[u, v, w] = getnc('u', 'v', 'wl', lims, lev, tim);
[lat, lon, depth] = getll('u', lims);
[latw, lonw, depthw] = getll('wl', lims);
u = u/100; v = v/100; w = w/1e6;
unan = find(isnan(u));
vnan = find(isnan(v));
u(isnan(u)) = 0;
v(isnan(v)) = 0;

[ntim, nlev, nlat, nlon] = size(u);

clear div
for i = 1:ntim;
  for j = 1:nlev;
    [div(i,j,:,:), lat2, lon2]...
         = sph_div1(squeeze(u(i,j,:,:)), squeeze(v(i,j,:,:)), lat, lon, 1, 1);
  end;
end;
div = shiftdim(div, 2);
div = [NaN * ones(1, nlon, ntim, nlev); div];
div = shiftdim(div, 1);
div = [div(nlon,:,:,:); div(1:(nlon-1),:,:,:)];
div = shiftdim(div, 1);
lon2 = [lon2(nlon)-360; lon2(1:(nlon-1))];

%  Integrate vertically

dz = diff([0; depthw]);
div = reshape(shiftdim(div, 1), nlev, nlat*nlon*ntim);
div = div .* (dz * ones(1, nlat*nlon*ntim));
div = cumsum(div);

div = shiftdim(reshape(div, nlev, nlat, nlon, ntim), 3);

%  Check the difference out

default_global; XAX = lonw; YAX = latw; FRAME = [0 360 -90 90];

figure(1); figure_orient;
tem1 = 1e6 .* squeeze(mean(w));
tem2 = 1e6 .* squeeze(mean(div));
tem2(tem2 == 0) = NaN;

lind = 2;
sp(1);
  gcont(squeeze(tem1(lind,:,:)), [-20:20]);
  dc2(squeeze(tem1(lind,:,:)));
sp(2);
  gcont(squeeze(tem2(lind,:,:)), [-20:20]);
  dc2(squeeze(tem2(lind,:,:)));

max(max(max(max(div - w))))/max(max(max(max(w))))
min(min(min(min(div - w))))/max(max(max(max(w))))

%  Now replace old values

tim = 46:146;  %  First check to see if the routine works
lev = 1:9;
lims = [-0.1 360 -90 90];

[u, v, w] = getnc('u', 'v', 'wl', lims, lev, tim);
[lat, lon, depth] = getll('u', lims);
[latw, lonw, depthw] = getll('wl', lims);
u = u/100; v = v/100; w = w/1e6;
unan = find(isnan(u));
vnan = find(isnan(v));
u(isnan(u)) = 0;
v(isnan(v)) = 0;

[ntim, nlev, nlat, nlon] = size(u);

clear div
for i = 1:ntim;
  for j = 1:nlev;
    [div(i,j,:,:), lat2, lon2]...
         = sph_div1(squeeze(u(i,j,:,:)), squeeze(v(i,j,:,:)), lat, lon, 1, 1);
  end;
end;
div = shiftdim(div, 2);
div = [zeros(1, nlon, ntim, nlev); div];
div = shiftdim(div, 1);
div = [div(nlon,:,:,:); div(1:(nlon-1),:,:,:)];
div = shiftdim(div, 1);
lon2 = [lon2(nlon)-360; lon2(1:(nlon-1))];

%  Integrate vertically

dz = diff([0; depthw]);
div = reshape(shiftdim(div, 1), nlev, nlat*nlon*ntim);
div = div .* (dz * ones(1, nlat*nlon*ntim));
div = cumsum(div);

div = shiftdim(reshape(div, nlev, nlat, nlon, ntim), 3);

cd /home/disk/hayes2/dvimont/csiro/data
nc = netcdf('wl_L2-10.nc', 'write');

  biff = nc{'wl'}(45,:,:,:);
  max(max(max(max(biff(biff <= 1e35)))))
  biff = nc{'wl'}(46,:,:,:);
  max(max(max(max(biff(biff <= 1e35)))))
  biff = nc{'wl'}(146,:,:,:);
  max(max(max(max(biff(biff <= 1e35)))))
  biff = nc{'wl'}(147,:,:,:);
  max(max(max(max(biff(biff <= 1e35)))))

     div = div * 1e6;
     mv = nc{'wl'}.missing_value;
     div(div == 0) = mv;
     nc{'wl'}(tim,lev,:,:) = div;

nc = close(nc);