cressman差值，程序、说明及例子
程序代码如下：
subroutine interp_cressman(lati,loni,data,nxi,nyi, !此处给定了nxi,nyi,nxo,nyo目的是给据网格上的i、j在主程序是方便计算经纬度的值，如经纬
* lato,lono,datao,nxo,nyo) !给定了 是否就不用再给出nxi,nyi,nxo,nyo呢？而只要给出待插值网格的循环次数呢？
!datao 是插值格点数据 !自己改写主程序和子程序的时候，可以使fvcom的经纬变成一维的L，就不必给出nxo,nyo了
c**** interpolates fields to another grid as specified
c**** by lat., lon. values.
c**** uses a Cressman interpolation with a specified
c**** search radius
c**** uses output from lookup.f to determine search points
parameter(pi=3.14159,num=500)
real lati(nxi,nyi),loni(nxi,nyi)
real lato(nxo,nyo),lono(nxo,nyo)
real data(nxi,nyi)
integer ilon(num),ilat(num)
real datao(nxo,nyo)
open(30,file='lookup.tab1',status='unknown')
open(31,file='lookup.tab2',status='unknown')
guess = 100. ! search radius
xmiss = -999.
scale = 1.
c************************************************************
c reinterpolate to NMC Octagonal grid using Cressman (应用cressman方法，逐步订正到NMC网格上)
c weights with a specified search radius
do 100 jj=1,nyo
do 100 ii=1,nxo
sum1=0.
sum2=0.
ilevs=0
c read in indices from lookup table
read(30,*) ilevs
c process indices
if(ilevs.ne.0.)then
read(31,*) (ilon(i),i=1,ilevs),(ilat(i),i=1,ilevs)
!ilon,ilat没有定义？ 有定义
do 90 kk=1,ilevs
i = ilon(kk)
j = ilat(kk)
c find distance between NMC point to be interpolated to
c and original point.
dim=(sin(pi*lati(i,j)/180.)*sin(pi*lato(ii,jj)/180.))
dam=(cos(pi*lati(i,j)/180.)*cos(pi*lato(ii,jj)/180.))
dam1=lono(ii,jj)-loni(i,j)
dam1=cos(pi*dam1/180.)
dist=dim+dam*dam1
dist=acos(dist)*110.949*(180./pi)
c dim=(sin(lati(i,j))*sin(lato(ii,jj)))
c dam=(cos(l
2021-11-23 19:37:51
1.78MB
1