修改6S Fortran77 代码，建立查找表-白红宇

修改6S Fortran77 代码，建立查找表

阅读量：4307 次

发布时间：2019-06-06

本文共 3392 字，大约阅读时间需要 11 分钟。

逐像元大气校正，常预先计算查找表（LUT，LookUp Tabel），6S大气辐射传输模式也可以用来计算LUT。但6S源程序输出信息多，且浮点数输出精度低，不利于提取关键信息生成LUT，本文描述了怎样修改6S源码以生成LUT。

首先确定LUT内容要素。

阅读MODIS M?D04 气溶胶产品生成算法文档（NASA相关网页），归纳了以下查找表要素：

1) 太阳天顶角观测天顶角太阳方位角观测方位角之差（相对方位角）散射角

2) 大气模式

3) 气溶胶模式

4) 550nm气溶胶光学厚度

5) 波段号

6) 大气透过率

7) atm. intrin. ref.（个人理解，这是大气程辐射换算后的反射率，用于校正计算）

8) total sca. （总散射，包括rayleigh散射和气溶胶散射，用于校正计算）

9) 表观反射率

10) 校正后的地表反射率

11) xa xb xc参数（物理意义不明，用于校正计算）

其次，阅读源码，明确LUT各要素在6S源码中的变量名。

6S大气校正计算源码（Excel验证中采用此公式）

rog=rapp/tgasm

rog=(rog-ainr(1,1)/tgasm)/sutott/sdtott

rog=rog/(1.+rog*sast)

xa=pi*sb/xmus/seb/tgasm/sutott/sdtott

xb=srotot/sutott/sdtott/tgasm

xc=sast

由计算源码确定需输出的变量。下面是输出LUT要素的Fortran77代码示例，建议放在源码main.f中stop一句之前。

write(*,*) asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55,

1 iwave,tgasm,ainr(1,1),sutott*sdtott,rapp,rog,xa,xb,xc

其中，asol是太阳天顶角，

phi0是太阳方位角，

avis是观测天顶角，

phiv是观测方位角，

adif是散射角，

phi是相对方位角，

idatm是大气模式号，

iaer是气溶胶模式号，

v是水平能见度，

taer55是550nm气溶胶光学厚度，

iwave表示波段号，

tgasm表示大气透过率，

ainr（1,1）是大气本身的反射率（姑且这么理解），

sutott*sdtott表示总散射，

rapp是表观反射率，

rog是校正后的地表反射率，

xa，xb，xc是校正计算参数。

Fortran77每行长度不能超过80个字符，续前行需在特定位置指明(示例中的iwave前的1即表示续行）。

该示例源码没有指定任何输出样式。浮点数会按默认样式输出，小数点后的位数比较多，精度较好。

挑选一个查找表文件用Excel验证后表明，excel公式计算值与6S输出值之间最大误差小于1E-7。说明方法是可行的。

再次，编译源码，编写Shell脚本：

编译环境：OpenSUSE操作系统 g95编译器，版本未知。

编译命令：g95 *.f -o sixs（在BRDF相关代码处可能有几个warning，本文不涉及BRDF，故暂不修改调试。在Windows下用f77编译，无warning，编译通过）

生成LUT的bash脚本getLUT.sh：

1:  function LUTCalc(){

2:  #{42,44,48,25,27,30}

3:  IBand=$1

4:  echo "Band ${IBand} is running"

5:  for Iref in {0.1,0.5}

6:  do

7:  fstat=${IBand}"_"${Iref}.res

8:  echo "asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55, iwave,tgasm,ainr,tott,rapp,rog,xa,xb,xc" >> ${fstat}

9:  for SolarZ in {0,15,30,45,75}

10:  do

11:    for SolarAz in {0,45,90,135}

12:    do

13:     for ViewZ in {0,15,30,45,75}

14:     do

15:      for Iaer in {1,2,3,5,6,7}

16:      do

17:       for Idatm in {1,2,3,4,5,6}

18:       do

19:        for IAod in {0.1,0.2,0.5,1.0};

20:        do

21:  #   Run sixs and output

22:  ./sixs <
     
      > ${fstat}

23:  0

24:  $SolarZ $SolarAz $ViewZ 0 3 15

25:  $Idatm

26:  $Iaer

27:  0

28:  $IAod

29:  -0

30:  -1000

31:  ${IBand}

32:  0

33:  0

34:  0

35:  0.0

36:  -$Iref

37:  EOF

38:         done

39:        done

40:       done

41:      done

42:     done

43:    done

44:    done

45:  }

46:  function getLUT()

47:  {

48:  echo "LUT is Calcing"

49:  for iwave in {42,44,48,25,27,30}

50:  {

51:    LUTCalc ${iwave}

52:  } &

53:  }

最后调用该脚本

>source getLUT.sh

>getLUT

最好晚上计算早晨看结果，如果CPU给力的话，几个小时后就可以得到结果。

下面是生成的LUT示例:

asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55, iwave,tgasm,ainr,tott,rapp,rog,xa,xb,xc

0.0000000 0.0000000 0.0000000 0.0000000 180.00000 0.0000000 1 1 63.664398 0.10000000 25 0.98984975 6.89081103E-02 0.80637234 0.10000000 3.87301818E-02 1.98730151E-03 8.71319696E-02 0.14777245

0.0000000 0.0000000 0.0000000 0.0000000 180.00000 0.0000000 1 1 26.739149 0.20000000 25 0.98984975 7.75793791E-02 0.76532620 0.10000000 2.94530466E-02 2.09388486E-03 0.10336593 0.16389242

0.0000000 0.0000000 0.0000000 0.0000000 180.00000 0.0000000 1 1 8.4940033 0.50000000 25 0.98984975 0.10173188 0.64870018 0.10000000 -2.69861287E-03 2.47033220E-03 0.15994170 0.20088956

0.0000000 0.0000000 0.0000000 0.0000000 180.00000 0.0000000 1 1 3.5674956 1.0000000 25 0.98984975 0.13688390 0.48083964 0.10000000 -7.89646432E-02 3.33272247E-03 0.29038188 0.24035060