计算Xt时维度也不一致

很抱歉，我再次犯了错误。计算Xt时的维度确实也存在问题。可以使用bsxfun函数或者直接使用矩阵相减的方式解决这个问题，具体实现如下：

使用bsxfun函数：

[a1, b1, ii] = ndgrid(1:jx, 1:jy, 1:zks);
JSX = a1.*gj - r0;
JSY = b1.*gj;
SLX = ZK(ii, 1);
SLY = ZK(ii, 2);
Xt = bsxfun(@minus, JSX, SLX);
rr = sqrt(sum(bsxfun(@power, [Xt, JSY - SLY], 2), 3));
[j, k] = ndgrid(1:nj, 1:nj);
z = (j - 1).*dz .* (j~=1 & j~=nj) + (j==1) + (j==nj).*(H-1);
a = (2*k-3).*dz*0.5 .* (k~=1 & k~=nj) + (k==1).*0 + (k==nj).*(H-dz/2);
b = (2*k-1).*dz*0.5 .* (k~=1 & k~=nj) + (k==1).*dz/2 + (k==nj).*H;

使用矩阵相减：

[a1, b1, ii] = ndgrid(1:jx, 1:jy, 1:zks);
JSX = a1.*gj - r0;
JSY = b1.*gj;
SLX = ZK(ii, 1);
SLY = ZK(ii, 2);
Xt = JSX - SLX(ii);
rr = sqrt(sum([Xt, JSY - SLY(ii)].^2, 2));
[j, k] = ndgrid(1:nj, 1:nj);
z = (j - 1).*dz .* (j~=1 & j~=nj) + (j==1) + (j==nj).*(H-1);
a = (2*k-3).*dz*0.5 .* (k~=1 & k~=nj) + (k==1).*0 + (k==nj).*(H-dz/2);
b = (2*k-1).*dz*0.5 .* (k~=1 & k~=nj) + (k==1).*dz/2 + (k==nj).*H;

这样就可以避免维度不一致的问题了。