粗化网格的变形速率计算

Calculate deformation rate of coarse-mesh

# 网格

# B-grid

CSIM uses a generalized orthogonal B-grid, where tracer quantities are located at the center of the grid cells, and velocities are located at the corners. The internal ice stress tensor takes four different values within a grid cell. Tracer quantities in the ice model include ice and snow area, volume, energy and temperature. The grid comprised of the center points of the grid cells is referred to as the “T grid”. The “U grid” is comprised of the ponts at the northeast corner of the corresponding cell on the T grid. Quantities that are defined on the U grid include ice and ocean dynamics variables.

计算应变速率所需要的速度是定义在 U grid 上的，也就是每个 cell 的东北角。论文中网格为$n_i$ 行$n_j$ 列，右端是高纬地区。

# 旋转

在截取了高纬地区 nj*5/7 : nj 列后，将下半部分旋转到右端拼接海洋。

% 更新维度为截取北极后

nj = size(uVel,2);

% 把下半数据旋转到右边

u = [ uVel([ni,1:ni/2],1:nj), -uVel(ni:-1:ni/2,nj-1:-1:1) ];

v = [ vVel([ni,1:ni/2],1:nj), -vVel(ni:-1:ni/2,nj-1:-1:1) ];

% 旋转后的网格去掉西边界和南边界

lat = [lat(1:ni/2,2:nj),   lat(ni-1:-1:ni/2,nj-1:-1:1) ];

lon = [lon(1:ni/2,2:nj),   lon(ni-1:-1:ni/2,nj-1:-1:1) ];

dx = [ dx(1:ni/2,1:nj), dx(ni:-1:ni/2+1,nj-1:-1:1) ];

dy = [ dy([ni,1:ni/2],2:nj), dy(ni:-1:ni/2,nj:-1:2) ];

% 更新维度为旋转后

ni = size(dx,1);

nj = size(dy,2);

新的网格:

# Strain rate 应变速率 /deformation rate 变形速率

根据论文中的定义，
$u_x=\frac{1}{A}\oint u \text{d}y$,
$u_y=-\frac{1}{A}\oint u \text{d}x$,
$v_x=\frac{1}{A}\oint v \text{d}y$,
$v_y=-\frac{1}{A}\oint v \text{d}x$,

剪切率 $\dot{\varepsilon}_{\textrm{shear}}=\sqrt{\left ( u_{x}-v_{y}\right )^2+\left ( u_{y}+v_{x}\right )^2}$,
散度率 $\dot{\varepsilon }_{\textrm{div}} = u_x+v_y$,
总应变率 $\dot{\varepsilon }_{\textrm{tot}}=\sqrt{\dot{\varepsilon }_{\textrm{shear}}^{2}+\dot{\varepsilon }_{\textrm{div}}^{2}}$。

在数值计算中，首先线性插值计算网格中点的速率（与$dx,dy$ 对应）

uMidy = (u(:,1:nj)+u(:,2:nj+1))/2;

vMidx = (v(1:ni,:)+v(2:ni+1,:))/2;

uMidx = (u(1:ni,:)+u(2:ni+1,:))/2;

vMidy = (v(:,1:nj)+v(:,2:nj+1))/2;

当网格粗化 times 倍时（以 3 为例），粗化后的网格单元即为上图中灰色的单元。

nx=fix(ni/times);  % 粗化后网格个数

ny=fix(nj/times);

ni=nx*times;

nj=ny*times;

每个大单元的速率偏导的计算公式为

$A=\sum_{i,j=1}^{3}dx_{i,j}dy_{i,j},$

$u_x \approx \frac{1}{A} \sum_{j=1}^{3} \left ( \left ( u_{midy_{4,j}} dy_{4,j} \right ) - \left ( u_{midy_{1,j}} dy_{1,j} \right ) \right ),$
$v_x \approx \frac{1}{A} \sum_{j=1}^{3} \left ( \left ( v_{midy_{4,j}} dy_{4,j} \right ) - \left ( v_{midy_{1,j}} dy_{1,j} \right ) \right ),$
$v_y \approx \frac{1}{A} \sum_{i=1}^{3} \left ( \left ( v_{midx_{i,4}} dx_{i,4} \right ) - \left ( v_{midx_{i,1}} dx_{i,1} \right ) \right ),$
$u_y \approx \frac{1}{A} \sum_{i=1}^{3} \left ( \left ( u_{midx_{i,4}} dx_{i,4} \right ) - \left ( u_{midx_{i,1}} dx_{i,1} \right ) \right ),$

对应代码：

% 计算每个大网格第 1：times 的相应数值

for i=1:times

    dxs(i,:,:)=dx(i:times:ni,1:times:nj);

    dxn(i,:,:) = dx(i:times:ni,1+times:times:nj+1);

    dye(i,:,:)=dy(1+times:times:ni+1,i:times:nj);

    dyw(i,:,:)=dy(times:times:ni,i:times:nj);

    ue(i,:,:) = uMidy(1+times:times:ni+1,i:times:nj);

    uw(i,:,:) = uMidy(times:times:ni,i:times:nj);

    ve(i,:,:) = vMidy(1+times:times:ni+1,i:times:nj);

    vw(i,:,:) = vMidy(times:times:ni,i:times:nj);

    vn(i,:,:) = vMidx(i:times:ni,1+times:times:nj+1);

    vs(i,:,:) = vMidx(i:times:ni,1:times:nj);

    un(i,:,:) = uMidx(i:times:ni,1+times:times:nj+1);

    us(i,:,:) = uMidx(i:times:ni,1:times:nj);

end

% 累加小网格的近似面积作为大网格的面积，经纬坐标取平均

for i=1:times

    for j=1:times

        area=area+((dx(i:times:ni,j:times:nj)+dx(i:times:ni,1+j:times:nj+1))/2).*...

            ((dy(i+1:times:ni+1,j:times:nj)+dy(i:times:ni,j:times:nj))/2);

        latMid=latMid+lat(i:times:ni,j:times:nj);

        lonMid=lonMid+lon(i:times:ni,j:times:nj);

    end

end

latMid=latMid/(times^2);

lonMid=lonMid/(times^2);

ux = ( ue .* dye - uw .* dyw );

vy = ( vn .* dxn - vs .* dxs );

uy = ( un .* dxn - us .* dxs );

vx = ( ve .* dye - vw .* dyw );

ux=squeeze(sum(ux,1)) ./ area;

vy=squeeze(sum(vy,1)) ./ area;

uy=squeeze(sum(uy,1)) ./ area;

vx=squeeze(sum(vx,1)) ./ area;

因此最后变形率的计算为：

div = ux + vy;

%shr = abs(ux-vy);

shr = sqrt((ux-vy).^2+(uy+vx).^2);

tot = sqrt(shr.^2+div.^2);