datareadnc.F

路径

LMDZ.MARS\libf\dynphy_lonlat\phymars\datareadnc.F

所属目录/模块

libf/dynphy_lonlat/phymars

文件定位

火星地表数据 NetCDF 读取与网格插值例程。datareadnc 从 surface.nc 文件读取地表粗糙度、反照率、热惯量、MOLA 地形和亚网格地形参数（hmons/summit/base），分别经水平插值或复合平均映射到 GCM 目标网格，输出给 newstart 用于初始化地表物理场。

定义的符号

依赖的模块

INCLUDE 头文件：

调用的关键例程

输入

输出

共享状态与副作用

• 读取 datafile_mod::datadir 模块变量，先赋默认值 "/u/lmdz/WWW/planets/mars/datadir"，再用 getin 从 callphys.def 覆盖。
• 打开并读取 surface.nc 文件（I/O 副作用）。
• 覆写 comconst_mod::pi：pi=2.*ASIN(1.)（行 111）。
• 无 common block 写入。

核心逻辑

1. 打开 NetCDF：设定 datadir 默认路径→getin 覆盖→NF_OPEN 打开 surface.nc；失败则打印指引信息并 CALL ABORT。
2. 读坐标：读取 latitude(jmdp1) 和 longitude(imd) 数组，转换为弧度标量盒坐标 rlonud(imdp1) 和 rlatvd(jmd)。
3. 主循环 k=0..4：依次读取 z0、albedo、thermal、zMOL（地形）、zMOL（亚网格统计）：
   • k=4（亚网格地形统计）：zMOL×1000 转 m → grid_noro1 计算亚网格统计量（zmea/zstd/zsig/zgam/zthe） → avg_horiz 计算 GCM 网格平均 → zavg。
   • k=0..3：数据填充到 zdataS（imdp1×jmdp1）→ interp_horiz 插值 → 周期性边界 → 按 k 保存到 z0（×0.01 转 m）/ alb / ith / phisinit。
4. 地形后处理：phisinit × 1000 × g 转换为地表位势（m²/s²）。
5. 扩展循环 k=5..7：读取 hmons、summit、base，使用 mvc_horiz 取最大值复合 → -999999 缺失值填 0。

伪代码

SUBROUTINE datareadnc(relief, phisinit, alb, ith, z0,
                      zmea, zstd, zsig, zgam, zthe,
                      hmons, summit, base, zavg)

  datadir ← "/u/lmdz/WWW/planets/mars/datadir"
  getin("datadir", datadir)
  NF_OPEN(datadir/surface.nc)
  READ latitude, longitude
  转换 → rlonud, rlatvd（弧度标量盒坐标）

  FOR k = 0 TO 4:
    string = ['z0', 'albedo', 'thermal', 'zMOL', 'zMOL']
    READ NetCDF variable string(k) → zdata
    IF k == 4 THEN  ! 亚网格地形
      zdata × 1000
      grid_noro1 → zmea/zstd/zsig/zgam/zthe
      avg_horiz → zavg
    ELSE
      zdata → zdataS（填充 imdp1 列）
      interp_horiz(zdataS → pfield)
      周期性边界补齐
      按 k 保存: z0(×0.01), alb, ith, phisinit
    ENDIF
  ENDFOR

  phisinit ← phisinit × 1000 × g

  FOR k = 5 TO 7:
    string = ['hmons', 'summit', 'base']
    READ NetCDF variable → zdata
    zdata → zdataS
    mvc_horiz（取最大值）→ pfield
    -999999 缺失值 → 0
    保存到 hmons/summit/base
  ENDFOR
END SUBROUTINE

参与的主题流程

写法特点

• Fixed-form Fortran（72 列限制、c 注释、&/$ 续行），但使用 ! 行注释和 :: 声明等 F90 扩展。
• NC_DOUBLE 条件编译：#ifdef NC_DOUBLE 选择 NF_GET_VAR_DOUBLE 或 NF_GET_VAR_REAL，适配不同 NetCDF 精度。
• 硬编码网格参数：imd=360, jmd=179, klatdat=180, ngridmxgdat=360 对应 surface.nc 的 360×180 全球 1° 网格。
• 缺失值处理：-999999. 作为 NetCDF 数据缺失标记，在扩展循环中填为 0。
• pi 覆写：行 111 pi=2.*ASIN(1.) 直接从反三角函数计算，不依赖外部赋值。

复现要点

• 需要 surface.nc 文件在 datadir 路径下可用（默认 /u/lmdz/WWW/planets/mars/datadir，可在 callphys.def 中覆盖）。
• 需要 dimensions.h 提供 iim, jjm 等 GCM 网格维度。
• relief='pla' 时地形场清零（平面模式），否则自动设为 'MOL'（MOLA 地形）。
• comgeom.h 提供 rlonu, rlatv 目标网格坐标。

待确认

• initial0 和 interp_horiz 所在文件/模块未在源码中显式 USE 或 EXTERNAL 声明，可能通过 include 或链接时接口引入。
• mvc_horiz_mod 的具体算法（最大值复合还是质量权重平均）需查看源码确认。

相关页面

• avg_horiz_mod — 水平空间平均插值，由本例程在 k=4 时调用
• grid_noro1 — 无旋转网格设置，由本例程在 k=4 时调用
• datafile_mod — datadir 模块变量来源
• dynphy_lonlat/index — 所属目录概览
• newstart — 调用方，在 start_archive 路径中读取 surface.nc 地表静态场