Cressman interpolation method

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(lati(i,j))*cos(lato(ii,jj))) c dam1=lono(ii,jj)-loni(i,j) c dam1=cos(dam1) c dist=dim+dam*dam1 c dist=acos(dist)*110.949 c if distance is .lt. 100 km then identify those c cases and find the cressman weights if(dist.le.guess .and. data(i,j).ne.xmiss)then zpiggy=(guess**2)-(dist**2) zpiggy1=(guess**2)+(dist**2) sum=zpiggy/zpiggy1 !(权重系数) c ilevs=ilevs+1 sum1=sum1+(sum*data(i,j)) !（权重系数乘以要素加和） sum2=sum2+sum !（权重系数加和） endif 90 continue endif c find interpolated value if(ilevs.ge.1) then datao(ii,jj)=(sum1/sum2)*scale !(订正系数) else write(*,*) 'Insufficient search radius ...exiting' datao(ii,jj)=xmiss endif 100 continue close(30) close(31) return end