1. 项目概述：为什么数组是Fortran的“王牌”？

如果你刚开始接触Fortran，可能会觉得这门语言有点“老古董”。但当你真正用它来处理科学计算、数值模拟或者大规模数据处理时，你会立刻发现它的魅力所在。而数组，就是Fortran这门“老将”手中最锋利、最趁手的武器。我最初接触Fortran是为了处理一个气候模型的海量网格数据，当时用其他语言写的脚本跑起来慢得像蜗牛，内存占用还高得吓人。后来切换到Fortran，仅仅是利用了其原生的、高效的数组操作，性能就提升了不止一个数量级。这让我深刻体会到，在Fortran的世界里，不懂数组，就等于没入门。

简单来说，数组就是一组具有相同类型的数据元素的集合，通过一个名字和下标来访问。这听起来和C、Python里的数组或列表差不多，对吧？但Fortran的数组远不止于此。它从语言设计之初，就是为了高效处理大型数值数组而生的。这种“基因”层面的优化，体现在内存布局、向量化操作、内置函数等方方面面。无论是处理一个简单的向量，还是一个庞大的四维时空数据场（比如温度、压力、速度在三维空间加一维时间上的分布），Fortran都能提供极其简洁和高效的语法支持。

所以，这篇教程的目标很明确：不仅要让你知道Fortran数组怎么声明、怎么用，更要让你理解它为什么这么设计，以及如何利用这些特性写出既快又稳的代码。无论你是要用Fortran做有限元分析（比如配合Abaqus）、计算流体力学，还是处理实验数据，数组都是你必须跨过去的第一道，也是最重要的一道坎。

2. 数组的核心概念与声明：从一维到多维

2.1 数组的声明与维度

在Fortran中，声明数组的核心是明确三要素： 数组名 、 数据类型 和 形状（维度与范围） 。最经典的声明方式是在变量名后加上用括号括起来的维度信息。

! 声明一个长度为10的一维实型数组
real :: temperature(10)

! 声明一个3行4列的二维整型数组
integer :: matrix(3, 4)

! 声明一个三维数组，例如表示一个空间体（x, y, z）上的物理量
real :: density(100, 100, 50)

这里需要注意下标范围的默认约定。在早期的Fortran（如Fortran 77）中，数组下标默认从 1 开始。这是Fortran的一个传统，也符合很多数学和物理公式的书写习惯（第一个元素是第1个）。在现代Fortran（Fortran 90及以后）中，你可以显式地指定下标范围，这带来了更大的灵活性。

! 显式声明下标从0开始，到9结束，共10个元素
real :: vector(0:9)

! 下标从-5到5，共11个元素，常用于对称区间处理
integer :: symmetric_index(-5:5)

! 二维数组，行从1到3，列从0到2
real :: grid(1:3, 0:2)

为什么下标从1开始是传统？ 这很大程度上源于其科学计算背景。在数学中，向量和矩阵的第一个元素通常标记为1（如a₁, M₁₁）。Fortran遵循这一惯例，使得代码公式和数学表达式之间的转换更为直观，减少了因下标偏移（如C语言的0起始）导致的“差一错误”（off-by-one error）。

2.2 动态数组：运行时决定大小

静态数组的大小在编译时就必须确定。但在实际科研中，我们经常需要处理大小在程序运行时才能确定的数据，比如从文件读取的未知行数数据集。这时就需要动态数组（Allocatable Arrays）。

声明动态数组使用 allocatable 属性，并在声明时用冒号 : 来标记未指定的维度。数组的实际大小通过 allocate 语句在运行时分配，使用完毕后应使用 deallocate 语句释放内存，这是一个好习惯。

program dynamic_array_demo
    implicit none
    real, allocatable :: data(:) ! 一维动态数组
    real, allocatable :: image(:, :) ! 二维动态数组
    integer :: n, m, i, stat

    ! 从用户输入或文件读取数组大小
    print *, '请输入一维数组的长度：'
    read *, n

    ! 分配内存，并检查是否成功
    allocate(data(n), stat=stat)
    if (stat /= 0) then
        print *, '内存分配失败！'
        stop
    end if

    ! 使用数组...
    data = 3.14

    ! 分配二维数组
    m = 5
    allocate(image(n, m))

    ! ... 使用 image 数组进行计算

    ! 计算结束后，释放内存
    deallocate(data, image)

end program dynamic_array_demo

实操心得：务必检查分配状态 。 allocate 语句的 stat 参数至关重要，特别是在处理可能非常大的数组时。内存不足是运行时常见错误，通过检查 stat ，程序可以优雅地失败并给出提示，而不是直接崩溃。养成 allocate 后必跟 if (stat/=0) 检查的习惯。

2.3 数组的整体操作与片段操作

这是Fortran数组语法中最优雅、最高效的部分之一。你不需要写循环就能对整个数组或数组的一部分进行操作。

整体赋值与运算：

real :: A(10, 10), B(10, 10), C(10, 10)

! 将数组A所有元素赋值为1.0
A = 1.0

! 数组对应元素相加，结果存入C。注意：形状必须完全一致！
C = A + B

! 数组与标量运算：所有元素乘以2.0
B = 2.0 * A

! 内置函数作用在整个数组上，返回一个同形状的数组
C = sin(A) ! 计算A中每个元素的正弦值

数组片段（Array Sections）： 你可以使用冒号 : 操作符来选取数组的一个连续子集，也就是片段。

integer :: vec(100)
integer :: mat(10, 10)

! 选取第5到第15个元素
vec(5:15) = 0

! 选取所有行（用冒号:表示），第3列
mat(:, 3) = 1 ! 将整个第三列赋值为1

! 选取第2到第5行，第4到第7列
mat(2:5, 4:7) = 9

! 使用步长，选取奇数索引的元素
vec(1:100:2) = -1 ! 第1,3,5,...,99个元素被赋值为-1

为什么这很强大？ 这种语法不仅让代码更简洁、更接近数学表达（比如 C = A + B ），更重要的是，它给了编译器极大的优化空间。编译器可以识别出这是对连续内存块的操作，从而更容易利用现代CPU的SIMD（单指令多数据流）指令进行向量化，大幅提升计算性能。你手写一个循环，编译器可能还需要分析才能决定是否向量化，而这种整体操作几乎是直接告诉编译器：“这里可以向量化！”

3. 数组在内存中的布局与高效访问

理解数组在内存中是如何存储的，是写出高性能Fortran代码的关键。Fortran使用 “列优先” 存储顺序。

对于一个二维数组 A(m, n) ，它在内存中是按列连续存放的：先存储第一列的所有行元素 A(1,1), A(2,1), ..., A(m,1) ，然后是第二列 A(1,2), A(2,2), ..., A(m,2) ，依此类推。

! 假设有数组 A(3, 2)
! 内存布局顺序是：
! A(1,1) -> A(2,1) -> A(3,1) -> A(1,2) -> A(2,2) -> A(3,2)

这与C/C++、Python (NumPy默认) 的“行优先”顺序正好相反。这个差异在混合编程（如用C调用Fortran库）时至关重要。

对性能的影响：访问局部性 现代CPU从内存中读取数据时，并不是一次只拿一个数，而是会一次性抓取一个“缓存行”（通常64字节）的数据到高速缓存中。如果程序访问的数据在内存中是连续的，那么缓存命中率就高，速度飞快；如果程序跳跃式访问，就会导致大量的“缓存未命中”，CPU不得不等待慢速的内存，性能急剧下降。

因此，在Fortran中遍历数组时，为了获得最佳性能， 内层循环应该是第一个下标（行索引） 。

! 高效的遍历方式（内层循环变行i）
real :: arr(1000, 1000)
integer :: i, j
do j = 1, 1000 ! 外层循环：列
    do i = 1, 1000 ! 内层循环：行
        arr(i, j) = arr(i, j) * 2.0
    end do
end do

! 低效的遍历方式（内层循环变列j）
do i = 1, 1000 ! 外层循环：行
    do j = 1, 1000 ! 内层循环：列
        arr(i, j) = arr(i, j) * 2.0 ! 每次访问都跳过了大量内存
    end do
end do

在低效的方式中，当 i=1 时，我们访问 arr(1,1) ，下一个访问 arr(1,2) 。但在内存中， arr(1,2) 离 arr(1,1) 很远（中间隔着第1列的其他999个元素），这破坏了访问的局部性。

注意 ：这是一个非常经典的性能陷阱。很多从C/C++转过来的程序员会不自觉地写出低效的循环顺序。记住口诀：“Fortran列优先，内层循环变行索引”。

4. 强大的内置数组函数

Fortran提供了一系列内置函数来处理数组，这些函数很多都能自动进行向量化并行计算，用起来非常方便。

4.1 形状查询与操作函数

• size(array [, dim]) : 返回数组的总元素数，或指定维度 dim 上的大小。

  integer :: A(5, 10)
  print *, size(A)     ! 输出 50
  print *, size(A, dim=1) ! 输出 5 (第一维大小)
  print *, size(A, dim=2) ! 输出 10 (第二维大小)
  

• shape(array) : 返回一个一维数组，包含各维度的大小。

  print *, shape(A) ! 输出 [5, 10]
  

• lbound(array [, dim]) / ubound(array [, dim]) : 返回数组在指定维度 dim 上的下标下界和上界。

  real :: B(0:4, -1:8)
  print *, lbound(B, dim=1), ubound(B, dim=1) ! 输出 0 4
  print *, lbound(B, dim=2), ubound(B, dim=2) ! 输出 -1 8
  

4.2 数学统计函数

这些函数通常可以带一个可选的 dim 参数。如果不指定 dim ，则对整个数组操作；如果指定了 dim ，则沿着该维度操作，结果数组的维度会减少一维。

• sum(array [, dim, mask]) : 求和。

  real :: vals(3) = [1.0, 2.0, 3.0]
  print *, sum(vals) ! 输出 6.0
  
  real :: mat(2,2) = reshape([1.,2.,3.,4.], [2,2])
  ! 矩阵： [1.0 3.0]
  !        [2.0 4.0]
  print *, sum(mat, dim=1) ! 按列求和，输出 [3.0, 7.0]
  print *, sum(mat, dim=2) ! 按行求和，输出 [4.0, 6.0]
  

• product(array [, dim]) : 求乘积。
• maxval(array [, dim, mask]) / minval(array [, dim, mask]) : 求最大值/最小值。

  print *, maxval(vals) ! 输出 3.0
  print *, minval(mat)  ! 输出 1.0
  

• maxloc(array [, mask]) / minloc(array [, mask]) : 返回最大值/最小值 所在的位置 （下标）。返回的是一个一维数组，长度等于原数组的秩（维度数）。

  integer :: pos(2)
  pos = maxloc(mat) ! 找出mat中最大值的位置
  print *, pos      ! 输出 [2, 2] (即第2行第2列，值为4.0)
  

• dot_product(vector_a, vector_b) : 计算两个一维数组（向量）的点积。

  real :: v1(3)=[1.,2.,3.], v2(3)=[4.,5.,6.]
  print *, dot_product(v1, v2) ! 输出 32.0 (1*4+2*5+3*6)
  

• matmul(matrix_a, matrix_b) : 矩阵乘法。这是Fortran进行线性代数计算的核心函数之一，高度优化。

  real :: A(2,3), B(3,2), C(2,2)
  A = reshape([1,2,3,4,5,6], [2,3])
  B = reshape([7,8,9,10,11,12], [3,2])
  C = matmul(A, B)
  ! C 将是 A 和 B 的矩阵乘积
  

实操心得：善用 mask 参数 。很多统计函数（如 sum , maxval ）支持 mask 参数，这是一个逻辑型数组，形状与原数组相同。只有对应 mask 为 .true. 的元素才会被纳入计算。这在处理带有无效值（如NaN）或需要条件统计的数据时非常有用。

real :: data(100)
logical :: valid(100)
! ... 给 data 赋值，并设置 valid 掩码（例如，标记大于0的数据为有效）
valid = (data > 0.0)
print *, '有效数据的平均值：', sum(data, mask=valid) / count(valid)

5. 数组的输入输出与格式化

将数组数据读写到文件或屏幕是常见操作。Fortran提供了灵活的输入输出方式。

5.1 隐式Do循环输出

这是一种简洁的输出数组所有元素的方式。

integer :: arr(5) = [1, 2, 3, 4, 5]

! 输出所有元素，用空格分隔
print *, arr
! 输出： 1 2 3 4 5

! 更可控的格式化输出，每行3个元素
write(*, '(3I5)') arr
! 输出：
!     1    2    3
!     4    5

(3I5) 是格式描述符，表示每行输出3个整数，每个整数占5个字符宽度。

5.2 读写文件

将数组写入文件或从文件读入同样方便。

program io_demo
    implicit none
    real :: data(100, 50)
    integer :: i, j, unit

    ! 生成一些示例数据
    do j = 1, 50
        do i = 1, 100
            data(i, j) = sin(0.1*i) * cos(0.05*j)
        end do
    end do

    ! 打开文件进行写入
    open(newunit=unit, file='output_data.dat', status='replace', action='write')
    ! 将整个二维数组写入文件，按列优先顺序
    write(unit, *) data
    close(unit)

    ! 从文件读回数据到另一个数组
    real, allocatable :: data_in(:, :)
    open(newunit=unit, file='output_data.dat', status='old', action='read')
    ! 动态分配数组大小需要先知道大小，这里假设已知
    allocate(data_in(100, 50))
    read(unit, *) data_in
    close(unit)
end program io_demo

注意 ：使用 write(unit, *) 这种列表导向格式输出时，数据是以空格分隔的纯文本。读取时， read(unit, *) 会按顺序读入足够的数据填满 data_in 。这种方式简单，但文件可能很大，且不保留维度信息。对于需要保存维度或更复杂结构的数据，可以考虑使用未格式化（二进制）文件或NetCDF、HDF5等科学数据格式。

5.3 格式化控制

对于需要精细控制输出格式的场景，格式描述符是关键。

real :: results(5) = [3.14159, 1.41421, 2.71828, 1.61803, 0.57721]

! 以表格形式输出，每个数占12个字符宽度，保留6位小数
write(*, '(5F12.6)') results
! 输出类似：
!    3.141590   1.414210   2.718280   1.618030   0.577210

! 每个数单独一行，带说明文字
do i = 1, 5
    write(*, '(A, I1, A, F10.5)') 'Result(', i, ') = ', results(i)
end do

6. 数组实战：一个简单的数值微分示例

让我们用一个完整的例子，把上面讲的知识点串起来。假设我们有一个由函数 f(x) = sin(x) 在区间 [0, 2π] 上离散采样得到的一维数组，我们要用中心差分法计算其导数的近似值。

中心差分公式： f'(x_i) ≈ (f(x_{i+1}) - f(x_{i-1})) / (2 * h) ，其中 h 是采样间距。

program numerical_derivative
    implicit none
    integer, parameter :: n = 100 ! 采样点数量
    real, parameter :: pi = 3.141592653589793
    real, parameter :: x_start = 0.0, x_end = 2.0 * pi
    real :: x(n), f(n), dfdx_approx(n)
    real :: h
    integer :: i

    ! 1. 生成等间距网格点和函数值
    h = (x_end - x_start) / real(n - 1)
    do i = 1, n
        x(i) = x_start + real(i-1) * h
        f(i) = sin(x(i))
    end do

    ! 2. 使用中心差分法计算内部点的导数
    ! 注意：边界点（i=1和i=n）无法用中心差分，这里简单处理
    dfdx_approx(1) = 0.0 ! 前向或后向差分，此处赋0
    dfdx_approx(n) = 0.0
    do i = 2, n-1
        dfdx_approx(i) = (f(i+1) - f(i-1)) / (2.0 * h)
    end do

    ! 3. 输出结果（前10个点）
    print *, '   x          sin(x)      approx cos(x)   exact cos(x)'
    print *, '--------------------------------------------------------'
    do i = 1, 10
        write(*, '(4F12.6)') x(i), f(i), dfdx_approx(i), cos(x(i))
    end do

    ! 4. 计算并输出整体误差（L2范数）
    ! 忽略边界点，计算内部点的误差均方根
    real :: error, rms_error
    error = 0.0
    do i = 2, n-1
        error = error + (dfdx_approx(i) - cos(x(i)))**2
    end do
    rms_error = sqrt(error / real(n-2))
    print *, ' '
    print *, 'RMS error (internal points): ', rms_error

end program numerical_derivative

这个例子涵盖了数组的声明、循环填充、基于下标的计算、格式化输出，以及简单的误差分析。你可以尝试修改函数 f(x) ，或者使用更复杂的差分格式。

7. 常见问题与避坑指南

在实际使用Fortran数组时，下面这些“坑”我几乎都踩过，希望你能避开。

7.1 数组越界（Index out of bounds）

这是最常见的运行时错误。访问了数组声明范围之外的下标。

integer :: arr(10)
arr(11) = 5 ! 编译可能通过，但运行时会崩溃或产生不可预知结果。

排查与预防 ：

1. 开启编译器边界检查 。在编译时加上调试选项，如gfortran的 -fcheck=bounds ，Intel Fortran的 /check:bounds 。这会在运行时捕获越界访问，并给出明确的错误信息。
2. 仔细检查循环的起始和结束条件。特别是当循环变量依赖于其他计算时。
3. 使用 lbound 和 ubound 函数来获取数组边界，而不是硬编码数字。

   do i = lbound(arr, 1), ubound(arr, 1)
       ! 安全操作
   end do
   

7.2 形状不匹配（Shape mismatch）

在进行数组整体赋值或运算时，等号两边的数组形状必须一致。

real :: A(5), B(10)
A = B ! 错误！形状 (5) 和 (10) 不匹配。

排查 ：编译器通常能捕获这类错误。确保进行操作的数组维度、各维度大小完全相同。使用 shape() 函数打印数组形状进行调试。

7.3 动态数组未分配或已释放

访问一个尚未分配（ allocate ）或已经释放（ deallocate ）的动态数组会导致严重错误。

real, allocatable :: temp(:)
temp(1) = 10.0 ! 错误！temp尚未分配内存。

allocate(temp(100))
! ... 使用 temp
deallocate(temp)
temp(50) = 20.0 ! 错误！temp已被释放。

预防 ：

1. 在使用前，用 allocated() 函数检查数组是否已分配。

   if (.not. allocated(temp)) then
       allocate(temp(100))
   end if
   

2. 释放数组后，最好将其显式置为空（ nullify 不适用于普通可分配数组，但可以将其指向一个空状态，或者更简单地，在释放后避免再次使用）。更好的做法是，在子程序或作用域结束时才释放数组。

7.4 性能陷阱：低效的循环与临时数组

除了前面提到的内存访问顺序问题，过度使用临时数组也会拖慢程序。

! 低效：创建了临时数组保存中间结果
real, allocatable :: tmp(:)
tmp = A + B
C = tmp * D

! 高效：合并操作，减少内存分配和拷贝
C = (A + B) * D

现代Fortran编译器很智能，像 (A+B)*D 这样的表达式，编译器通常会生成优化的代码，避免生成完整的临时数组 A+B ，而是采用融合循环的方式计算。但如果你显式地分配了 tmp ，就强制了中间结果的存储。

7.5 与C语言交互时的顺序问题

当你用C语言调用Fortran子程序，或者反过来时，必须特别注意数组的存储顺序。

假设Fortran有一个子程序：

! Fortran side: sub.f90
subroutine process_matrix(mat, m, n)
    integer, intent(in) :: m, n
    real, intent(inout) :: mat(m, n) ! 列优先
    ! ... 处理矩阵
end subroutine process_matrix

在C语言中调用时，传入的数组必须是 转置后 的，或者C代码需要按照列优先的顺序来理解和填充数据。

/* C side: main.c */
extern void process_matrix_(float *mat, int *m, int *n); /* 注意名称可能加了下划线 */

int main() {
    int m=3, n=2;
    float mat[3][2]; /* 在C中是行优先: mat[行][列] */
    /* 填充数据... */
    /* 但Fortran期望的是列优先，所以这里直接传递可能导致数据解释错误 */
    process_matrix_(&mat[0][0], &m, &n);
    return 0;
}

解决方案 ：要么在C端将数据按列优先存储，要么在接口层进行转置。对于复杂的交互，建议使用ISO_C_BINDING模块，它提供了更安全、标准的互操作方式。

数组是Fortran的灵魂，掌握它，你就掌握了用Fortran进行高效数值计算的钥匙。从简单的声明到高效的内存访问，从强大的内置函数到实际的应用案例，每一步都需要理解和练习。刚开始可能会觉得有些概念（如下标顺序、片段操作）比较别扭，但一旦习惯，你就会爱上这种直接、高效的表达方式。我个人的体会是，在处理真正的多维网格数据时，Fortran数组的简洁性和性能是其他语言难以比拟的。多写，多试，多调优，尤其是多用编译器提供的诊断和性能分析工具，你会进步飞快。