程序的调用过程

类似下面链接中介绍的调用步骤，在LMDZ.MARS/libf/phymars/physiq_mod.F调用physiq子程序，会进一步调用辐射传输计算子程序callradite（callradite_mod.F），在其中进一步调用子程序aeropacity（aeropacity_mod.F）。

数据的存储

physiq子程序的输入为基本参数以及各物理量，沙尘信息存在于tracers变量pq中。整个子程序的输出为基本物理量的倾向，沙尘倾向存于变量pdq中。
callradite子程序的输入为pq（在physiq子程序中同变量名），主要输出有tauscaling, flux…, tau, aerosol等（同变量名）
callradite子程序中，除了用于输出外，tauscaling和tau只在调用子程序updatereffrad和aeropacity中出现，前者似乎是用于计算有效半径。aerosol则在aeropacity中出现（完全由aeropacity计算出的），其后用于计算zaerosol，而zaerosol仅用于计算短波与长波辐射，也不用于输出。

<aside> 💡 tauscaling和tau完全由aeropacity计算出的，但是在调用updatereffrad时，采用的是上一次的计算值或初始文件读入的预设值。

</aside>

如果在def文件中设置了rdstorm变量（rocket duststorm）的值为true（默认为false）

Subroutine aeropacity

在LMDZ.MARS/libf/phymars/aeropacity_mod.F文件中包含子程序aeropacity。

程序的第一步是根据读入的不同iaervar值，读取相应的预设沙尘文件，将结果存入tau_pref_scenario中。（原始文件在309行附近）。具体方法是调用read_dust_scenario子程序。

#  =1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario,
#  =4 Mars Year 24 from TES assimilation (old version of MY24; dust_tes.nc file)
#  =24 Mars Year 24 from TES assimilation (ie: MCD reference case)
#  =25 Mars Year 25 from TES assimilation (ie: a year with a global dust storm)
#  =26 Mars Year 26 from TES assimilation

aeropacity程序的第二步是根据各个aerosol的名称，选择相应的方式计算出aerosol的值（根据传入的pq等变量）
aeropacity程序的第三步是对沙尘这种aerosol进行进一步处理。首先，会根据reff_driven_IRtoVIS_scenario变量，来决定是否重新读取一下预设沙尘文件（原文件703行附近）。对于默认设置，该变量值为false，故这一步会跳过。
然后，调用compute_dustscaling子程序，结合之前读入的tau_pref_scenario数据，给出修正系数（tauscaling或者dust_rad_adjust，目前用前者）和修正后的气溶胶数组aerosol。(过程见下一节）
接下去计算透明度tau_pref_gcm，以传出供其他程序使用。如果freedust的值为false（该值从def文件中读入，但是自带的def文件中未给出，故一般为默认值false），那么tau_pref_gcm就直接等于读入的tau_pref_scenario。如果freedust的值为true，则根据修正后的气溶胶数组aerosol计算。
aeropacity程序的第四步是通过对aerosol垂直求和给出变量tau，用于输出。
如果rdstorm或topflows的值为true（默认均为false），aeropacity程序最后还会对aerosol作一些调整。

Subroutine compute_dustscaling

在LMDZ.MARS/libf/phymars/dust_scaling_mod.F90文件中，包含子程序compute_dustscaling。根据dustscaling_mode变量的值（该值从def文件中读入，但是自带的def文件中未给出，故一般为默认值1），下面有三种逻辑。
若dustscaling_mode为0，则tauscaling直接为1，即不进行任何scale操作。

若dustscaling_mode为1，则先根据传入aerosol的变量计算每个水平网格的沙尘总量，存入taudust中。然后，根据如下公式计算tauscaling。式中odpref为常数610，表示沙尘光深选取的参考气压。pplev(:,1)为第一层inter-layer的气压。最后，对所有的沙尘种类，乘以上面计算出的tauscaling，存入aerosol数组。

taudust(:) = 0
      do iaer=1,naerdust ! loop on all dust aerosols
        do l=1,nlayer
          do ig=1,ngrid
              taudust(ig)=taudust(ig)+aerosol(ig,l,iaerdust(iaer))
          enddo
        enddo
      enddo
--------
tauscaling(:)=tau_pref_scenario(:)*pplev(:,1)/odpref/taudust(:)
--------
do iaer=1,naerdust
        do l=1,nlayer
          do ig=1,ngrid
            aerosol(ig,l,iaerdust(iaer)) = max(1E-20, &
                      aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
          enddo
        enddo
      enddo

若dustscaling_mode的值为2，则在计算taudust时，只采用一种沙尘（程序描述为background dust）。然后判断当前火星日是否已经计算过dust_rad_adjust，若没有，则调用子程序compute_dust_rad_adjust进行计算。最后，对所有的沙尘种类，乘以dust_rad_adjust，存入aerosol数组。
在子程序compute_dust_rad_adjust中，dust_rad_adjust的计算方式和上面tauscaling的公式类似，但是进一步根据当地的时间，计算相邻两天tauscaling的加权平均。

程序的调用过程

Subroutine aeropacity

Subroutine compute_dustscaling

注