SUNTANS模型学习（2）——用MATLAB工具提取Time Accuracy算例的计算数据

算例回顾（Time Accuracy算例）

该算例计算区域为垂向的二维立面。区域的水平长度为100m，深度为100m。
水平方向上，该计算域由一百个等边三角形网格组成；水平向网格示意图如下：

在垂直方向上，由50个相等的垂直网格组成；垂向网格示意图如下：

模拟设置与参数详见之前的博文SUNTANS模型学习（1）——运行Time Accuracy算例。

数据提取MATLAB工具

在SUNTANS模型的源代码包中，有一名为"mfiles"的文件夹，我们打开getcdata.m，看到如下代码：

%
% GETCDATA
% DATA = GETCDATA(FID,ARRAYSIZE,STEP,PRECISION);
%   extracts binary data from the specified file descriptor.  As
%   an example, if a file contains double precision arrays of length
%   N=10, the n=10th array is extracted with:
%   
%   data = getcdata(fid,N,n,'float64');
%
% Oliver Fringer
% Stanford University
% 27 July 05
%
function data = getcdata(fid, arraysize, step, precision)

  tic;
  if(fid>0),
    if(precision=='float32'),
      space = 4*arraysize;                    
    elseif(precision=='float64')
      space = 8*arraysize;
    end

    frewind(fid);
    fseek(fid,(step-1)*space,0);
    data = fread(fid,arraysize,precision);
  else
    fprintf('Error in function: getdata...undefined fid\n');
    data = zeros(arraysize,1);
  end
  
  fprintf('Read data at a rate of %.2f Mb/sec\n',space/1e6/toc);

我们提取数据时所用的就是上述getcdata函数。它的输入参数包括：

1. fid：数据文件（.dat）的句柄；
2. arraysize：每一组数据的个数；
3. step：时间步；
4. precision：精度，可选‘float32’或float64’；

示例：提取Time Accuracy算例数据

首先利用 addpath() 函数去调用 getcdata.m 函数：

PATH = 'D:\suntans-master\mfiles'; % PATH 为 getcdata.m 所在文件夹的绝对路径
addpath(PATH);

之后，选择MATLAB的工作区域为数据文件所在的文件夹；执行下列命令：

filename = 's.dat';	           
% 此处以 s.dat （盐度场数据文件）为例；盐度值储存于网格中心；
fid = fopen('s.dat','r+');
Ncell = 50;	Nkmax = 50;
% 由之前的算例设置可知，计算域每层共有50个网格，垂向共50层；
T_step = 2;
% 此处提取第2个时间步时盐度场数据；

data = getcdata(fid, Ncell*Nkmax, T_step, 'float64');
% 提取数据得到data，它的数据总个数会等于 总时间步数、每层的网格数和垂向层数之积；

最后将提取的一列数据加以整理：

salinity = zeros(Nkmax,Ncell);
for i = 1:Nkmax
p_start = Ncell*(i-1) + 1;
p_end = Ncell*i;
salinity(i,:) = data(p_start: p_end); 
end
% 得到的 salinity 矩阵即 T_step 时间步下的盐度场
% salinity 的每一行代表每一层的盐度；而每一列表示对应网格处垂向盐度分布（索引值越大，所表示的网格越靠近底部）；

最后画出垂直方向上的盐度分布图：

 pcolor(flip(salinity)),shading flat, colorbar, colormap jet