上周什么都没有做，就把那个micromet插值弄了一下，现在把它总结下方便后面有需要的孩纸，总结得不到位的还请原谅！

先来简单的介绍下micromet这个东东吧！MicroMet 是专为插值大气驱动而设计的半物理过程模型[126]，插值方法考虑到高程的变化，使用数学方法、经验公式和物理公式结合的办法对大气驱动进行插值，可插值八个驱动气象要素，包括气温、相对湿度、风速、风向、向下长波、短波辐射、气压和降水，另外模型内部还包含一个在时间序列上插值补充缺省值的计算方法。具体对各驱动变量插值时，MicroMet 第一步首先采用纯数学的距离加权插值方法将源数据插值到网格，并且将输入的源数据的高程信息也按距离加权插值出一个虚拟的高程水准面，第二步再对逐个网格参考网格真实高程与虚拟高程水准面间的高程差单独调整插值结果：插值气温时用高程差与气温直减率来计算；相对湿度则是先结合气温换算出露点温度，设定露点直减率，插值出露点，最终反算出比湿；插值风速时则是先把风速拆成南北向东西向两个分量，分别插值再合二为一，后用风速、坡度、曲率半径来作调整；插值向下长波、短波辐射时则是结合云反射率、温度和高程来调整；气压按高程换算。MicroMet插值降雨时采用了非线性的降雨—高程方程[127]，其中涉及一个重要的参数降水调整系数。在制备数据的过程中发现，降水调整系数选取的是否合适严重影响降水驱动插值结果。MicroMet 插值降水参考的公式为：

上面的内容摘自我师姐侯玉婷的毕业论文，如有引用，请你注明，哈哈！文章题目：黑河流域上游高分辨率降水驱动析及 Noah LSM的径流响应模拟！

看了上面的东西你是不是还是不是很懂这个呀，当然的，因为前面只是一个简单的说明，具体的文献请你参考：A Meteorological Distribution System for High-Resolution Terrestrial Modeling，这篇文章有详细的micromet插值原理，建议你在插值前细细的读读这篇文章。那么下面我就具有来说说这个插值的过程，我插的是降雨，因此下面主要针对降雨插值的介绍。

micromet是开源的，因此我也会毫不吝啬的发给你，在这里面的说明都被我翻译成中文了，希望没有翻译错，如果需要，请你邮件：zhanglingky@126.com，并说明你的目的，谢谢！

我的实际观测数据是站点降雨数据，由于我研究区域站点较少，不能很好的作为水文模型的驱动数据。因此非常有必要对驱动站点数据进行插值处理。下面看看我的的研究区域：

红色的站点是我的实际站点，而蓝色的站点是我需要降雨的站点，但是实际上这些站点没有数据，因此通过插值得到格网的降雨实际序列数据，然后用站点所在的格网数据代表站点数据。

    首先要准备的数据包括3个，DEM、植被数据（这个数据插降雨是不需要的，但是模型要求任然要给出，因此只需要对DEM重分类为1即可作为fake_veg）。还有就是站点数据。下面说数据准备。

一、站点数据的组织都有很好的说明，值得注意的是需要把这些站点数据综合起来，这个在文件中有一个专门的fortan程序来实现，名字为：mk_micromet_input_file.f 在文件夹met里面，顺便说一下由于micromet是用fortran写的，因此你大电脑必须安装fotran编译器，而且micromet 程序的调用需要使用到gfortran，这个也是你需要安装的。

二、Dem数据的大小决定你插值分辨率的大小，而且你的植被（假的）需要和dem分辨率和投影系统一样。这样才才能确保程序正常运行。数据准备好转化为ASC格式。放到文件夹topo_veg

三、前面两个工作准备了，现在就是模型的配置的问题，这个比较重要，用ultraedit打开snowmodel.par这个文件，里面有配置的说明，现在我说几个需要注意的地方，

! Number of x and y cells in the computational grid.
     nx = 50
     ny = 35

这个nx和ny是你ASC文件即dem文件的行列数，nx为列数，ny为行数。

! deltax = grid increment in x direction.  Meters.
! deltay = grid increment in y direction.  Meters.
     deltax = 5000.0
     deltay = 5000.0

! Location (like UTM, in meters) value of lower-left grid point.
!   xmn = value of x coordinate in center of lower left grid cell.
!     Meters.
!   ymn = value of y coordinate in center of lower left grid cell.
!     Meters.
     xmn =  -571811.53424883
     ymn = 4043727.467686

deltax为分辨率；xmn和ymn为asc文件左下角坐标。

! Model time step, dt.  Should be the same increment as in the input
!   data file.  In seconds.
!   One day.
    dt = 86400.0
!   Six hours.
!     dt = 21600.0
!   Three hours.
!     dt = 10800.0
!   One hour.
!      dt = 3600.0

这个确定你的插值的时间步长。我插的是天为步长的。

! Number of model iterations defines how many times to process.
!      max_iter = 7656
      max_iter = 7670

这个确定你插值的长短，迭代一次，为一次步长。

! Define whether the precipitation adjustment factor, with units of
!   km^-1 (kind of a precipitation lapse rate, used to adjust the
!   precipitation for locations above and below the precipitation
!   observing station(s)), is to use the default monthly lapse
!   rates (= 0) or user supplied monthly lapse rates (= 1).  To use
!   user supplied lapse rates, you have to edit the user lapse rate
!   data array in micromet_code.f (subroutine get_lapse_rates).
     iprecip_lapse_rate_user_flag = 1

这个参数设为1表示，用自己定义的降雨调节因子，这个因子的算法我后面讲。

设置好，下面就是模拟，在dos下打开代码所在的文件夹，输入built，完成配置；然后返回到上一级，输入snowmodel.exe模型就开始运行了。

运行结果放在output里面，这个结果存放格式为.gdat，

DSET  ^micromet.gdat
TITLE SnowModel
UNDEF -9999.0
XDEF   31 LINEAR 0 0.2
YDEF   31 LINEAR 0 0.2
ZDEF    1 LEVELS 1
TDEF  319 LINEAR 01Z01oct2002 1dy
VARS    9
ta          1  0 air temperature (deg C)
rh          1  0 relative humidity (%)
u           1  0 meridional wind component (m/s)
v           1  0 zonal wind component (m/s)
wspd        1  0 wind speed (m/s)
wdir        1  0 wind direction (0-360, true N)
qsi         1  0 incoming solar radiation reaching the surface (W/m2)
qli         1  0 incoming longwave radiation reaching the surface (W/m2)
prec        1  0 precipitation (m/time_step)
ENDVARS

一般的软件打不开，因此我师姐曹艳萍找人写了个代码把它读出来，很简单，你自己看看

   implicit none
integer nx,ny,nv,nt,var,xnum,ynum
parameter(nx=50,ny=35,nv=9,nt=7670,var=9,xnum=31,ynum=21)
character*16 str(9)
data str /'ta','rh','u','v','wspd','wdir','qsi','qli','prec'/
!real air(nx,ny,nv,nt)
real air(nx, ny, nt)
real tmp(nx, ny)
   integer i,j,v,t
open(1,file='F:micrometmicrometmicrometoutputsoutput.gdat',form='binary',status='old')
do t=1,nt
do v=1,nv
 if(v==var)then
  !read(1) ((air(i,j,v,t),i=1,nx),j=1,ny)
  read(1) ((air(i,j,t),i=1,nx),j=1,ny)
 else
  read(1) ((tmp(i,j),i=1,nx),j=1,ny)
 end if
end do
print*, 't=',t
end do
print*,' read date completed'
close(2)
close(1)

5 FORMAT(50f12.6)
100 format(f12.6)

   open(2,file='F:micrometmicrometmicrometoutputs'//trim(str(var))//'.txt',status='replace')
do t=1,nt
!write(2,5) ((air(i,j,var,t),i=1,nx),j=1,ny)
write(2,5) ((air(i,j,t),i=1,nx),j=1,ny)
end do
   close(2)

open(3,file='F:micrometmicrometmicrometoutputs/'//trim(str(var))//'_singlepoint.dat',form='formatted',status='replace')
do t=1,nt
!write(2,100) air(xnum,ynum,var,t)
write(3,100) air(xnum,ynum,t)
end do

end

红色部分是你需要修改的，nx,ny就不说了对应5 FORMAT(50f12.6) 50修改为你的nx。var表示你提取的变量，我提前降雨所以是9，因此后面2个xnum=31,ynum=21表示你提前的某一个格网的数据，31列和21行的数据。

值得注意的是，写出来的降雨上下和arcgis实际上下颠倒，而arcgis的行和实际的xy坐标中的行也是颠倒的，因此我们再提取某个位置的数据，就不需要对ynum进行了转化了，也就是说这里面的多少行多少列就是你实际arcgi中的多少行和多少列。

但是我们再提取面上的数据时候，就需要把它颠倒过来，这样才能得到实际的数据。我的某一年插值结果如下：

最后我说说降雨调节因子的算法。

任意2站点组合。比如你有6个站点，那么组合共6*5/2共15个方程，且每个方程都是超定方程，而所有组合合起来就是超定之超定方程，实际解放很简单，matlab左除的方法能够解决这样的超定方程，求得最小二乘解