修改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：

   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 <<EOF | tail -n 1 >> ${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   
•        ……

转载于:https://www.cnblogs.com/wishmo/p/3429484.html