修改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???
• ?????? ……