用 Matlab 计算 IGRAv2 探测站的 Tm 和 PWV

1. 探空站IGRAv2数据

探空站的Tm常作为真值，去检验Tm线性公式或者ERA5 Tm等的精度 。

探空站PWV常作为真值，去检验GNSS PWV等的精度

2. Tm 的计算方法

Tm 的计算方法有两种在前面的文章有讲，这里用
使用水汽压和温度计算Tm。

ei和 Ti 表示第 i 层大气的水汽压和温度,Zi是第 i 层的厚度。这种方法使用探空站数据计算得到的 Tm 一般为真值与其他数据进行比较。

2.1. 水汽压（WVP）计算公式

根据饱和水汽压(es)和相对湿度(RH)求得，公式如下

2.2. es 为饱和水汽压（hPa），用下式求得

式中，T 为温度（K）。es 采用 ECMWF IFS 报告（IFS Documentation CY31R1 Part II ）给出的模型，对水的不同状态做了区别 (ECMWF, 2007) ：  
   （1）温度大于 0℃， R2 = 611.21 hPa，R3 = 17.502 K 和 R4 = 32.19 K；  
   （2）温度小于-23℃，R2 = 611.21 hPa，R3 = 22.587 K 和 R4 = -0.7 K；  

    （3）温度介于-23℃ 和 0℃ 之间，则用下式计算：

                                         式中，T0 = 273.16 K，Ti = 250.16 K。

3. PWV 的计算方法

3.1. 计算比湿度q

e是水气压，计算方法在上面得Tm计算里面有。

3.2. 重力加速度g考虑了纬度和高程变化。

4.Matlab 获取代码关注公众号WZZHHH，或者咸鱼关注：WZZHHH123，部分代码展示，总共337行代码：
 

% 读取IGRA的文件
% 基于matlab2022写的
% 微信公众号和****：WZZHHH
% IGRA源数据下载地址https://www1.ncdc.noaa.gov/pub/data/igra/derived/derived-por/
% 里面有很多站点已经弃用，并且该网址包括建站以后的所有数据集,
% 因此我们需要提取对应年份的数据,然后计算得到探空站得PWV

% 提取数据存到"v2data"，按照时间进行排列，最后存到mat文件中
% 每个站点的每年存一个mat文件，名称为：站点名字-年份.mat
% --------------数据的含义如下所示---------------
%    V2data.pwv：大气可降水量;
%    V2data.lon：经度;
%    V2data.lat：纬度;
%    V2data.Tm ：加权平均温度;
%    V2data.ID ：站点名称;
%    V2data.pwv有6列：第1列是pwv，第2-6列是时间和时间数字

% ------------------------------------------
clc,clear;

%% -----------------------------需要修改的参数--------------------------
% 列 参数           文本位置           含义
%  1  PRESS           1-  7   Integer 压力(Pa或mb*100)
%  2  REPGPH          9- 15   Integer 报告的位势高度(米)。这个值在重要级别上是不可用。
%  3  CALCGPH        17- 23   Integer 计算出的位势高度(米)
%  4  TEMP           25- 31   Integer 报告的温度(K*10)。
%  5  TEMPGRAD       33- 39   Integer 当前高度与下一高度之间的温度梯度[(K/km)*10，如果温度随高度增加，则为正]。
%  6  PTEMP          41- 47   Integer 位势温度（K*10）
%  7  PTEMPGRAD      49- 55   Integer 当前高度和下一高度之间的位势温度梯度[（K/km）*10，如果潜在温度随高度增加，则为正]。
%  8  VTEMP          57- 63   Integer 虚温(K*10)
%  9  VPTEMP         65- 71   Integer 虚势温度(K*10)
%  10 VAPPRESS       73- 79   Integer 同一水平的温度、压力和露点下降计算得出的水气压
%  11 SATVAP         81- 87   Integer 饱和水汽压
%  12 REPRH          89- 95   Integer 为原始探测时报告的相对湿度(百分比*10)
%  13 CALCRH         97-103   Integer 为相对湿度(百分数* 10)，由蒸汽压、饱和蒸汽压和同一水平面的压力计算得到。
%  14 RHGRAD        105-111   Integer 当前高度和下一高度之间的相对湿度梯度[(%/km)*10，如果相对湿度随高度增加，则为正]。
%  15 UWND          113-119   Integer 纬向风分量[(m/s)*10]，由报告的风速和风向计算得到。
%  16 UWDGRAD       121-127   Integer 当前高度和下一高度之间纬向风的垂直梯度[(m/s/km)*10，纬向风随高度增大而增大时为正]。
%  17 VWND          129-135   Integer 经向风分量[(m/s)*10]，由报道的风速和风向计算得到。
%  18 VWNDGRAD      137-143   Integer 当前高度与下一高度之间经向风的垂直梯度[(m/s/km)*10，经向风随高度增大而增大则为正]。
%  19 N             145-151   Integer 折射率(无单位)

% igra_num：需要提取数据所在的列
% 计算PWV所需要得数据为[温度 气压 计算位势 计算相对湿度]
% 改对应第4 1 3 13 列；
igra_num = [4 1 3 13];

% IGRA路径
path_igra = 'D:\DATA\IGRA2\drvd\';

% 数据存放路径
save_igra = 'D:\DATA\IGRAv2PWV\';

% IGRA官网下载的站点清单（igra2-station-list）存放路径
% 清单下载地址：https://www1.ncdc.noaa.gov/pub/data/igra/igra2-station-list.txt
list_path = 'D:\paper_write\paper_code\2\IGRA\igra2-station-list.txt';

% IGRA数据包含了建站以来数据，挑选所需年份
year = [2020 2021]; % 以2020-2021年为例