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

 

逐像元大气校正，常预先计算查找表（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：

function LUTCalc(){
#{42,44,48,25,27,30}
IBand=$1
echo "Band ${IBand} is running"
for Iref in {0.1,0.5}
do
fstat=${IBand}"_"${Iref}.res
echo "asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55, iwave,tgasm,ainr,tott,rapp,rog,xa,xb,xc" >> ${fstat}
for SolarZ in {0,15,30,45,75}
do
  for SolarAz in {0,45,90,135}
  do
   for ViewZ in {0,15,30,45,75}
   do
    for Iaer in {1,2,3,5,6,7}
    do
     for Idatm in {1,2,3,4,5,6}
     do
      for IAod in {0.1,0.2,0.5,1.0};
      do
#   Run sixs and output
./sixs <<EOF | tail -n 1 >> ${fstat}
0
$SolarZ $SolarAz $ViewZ 0 3 15
$Idatm
$Iaer
0
$IAod
-0
-1000
${IBand}
0
0
0
0.0
-$Iref
EOF
       done
      done
     done
    done
   done
  done
  done
}
function getLUT()
{
echo "LUT is Calcing"
for iwave in {42,44,48,25,27,30}
{
  LUTCalc ${iwave}
} &
}

• 最后调用该脚本

>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   
•        ……