EC细网格风场格点数据插值为站点数据问题

hushunqi215

主程序：
！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！


!   双线性插值，格点到站点，考虑权重  !
!      modifioed from    mofangbao(清风) by djj!   

!===================================================!
program bilineargird2stn
integer,parameter::oldx=961,oldy=561    !格点数据xy向网格点数
integer,parameter::newxy=9             !站点数
integer,parameter::time=1               !时间层次
real,parameter::oldLon_start=10.0,oldLat_start=10.0,oldInterpx=0.125,oldInterpy=0.125  !格点数据起始经纬度和格距
real oldVal(oldx,oldy),newVal(newxy),newPoint(newxy,2)  !分别是：格点数据，站点数据，站点数据经纬度，newVal的点顺序和nowPoint中的对应
character*8 ::stid
integer::i,j
logical isFirst   !判断是否第一个时次
isFirst=.true.          !后面无需手动修改此值，程序会自动修改


!!读取站点经纬度
open(27,file='E:\yh\station.txt')
do i=1,newxy
        read(27,*) stid,newPoint(i,2),newPoint(i,1)  !先经度后纬度
enddo
close(27)


!!读取网格点数据
open(29,file='E:\yh\18052220.120',form='binary')
!read(29,*)
!read(29,*)
!read(29,*)
read(29,*)  ((oldVal(i,j),i=1,oldx),j=1,oldy)
close(29)



!!调用插值程序,这里只哟一个时次，因此没有循环
call getNewVal(oldVal,newVal,oldLon_start,oldLat_start,oldInterpx,oldInterpy,newPoint,isFirst)
open(25,file='E:\yh\grid2stn.txt')
do i=1,newxy
write(25,*) newVal(i)
enddo
end




subroutine getNewVal(oldVal,newVal,oldLon_start,oldLat_start,oldInterpx,oldInterpy,newPoint,isFirst)
integer,parameter::oldx=961,oldy=561,newxy=9
real::oldLon_start,oldLat_start
real,parameter::undef=-999.0  
real::oldInterpx,oldInterpy      
integer pos(newxy,4),nstart,nend  !weight 1,2是左侧点和右侧点的分别权重，3,4是下面点和上面点分别的权重
real tmpVal_up,tmpVal_down !本次纬向下方和纬向上方的线性插值结果以及权重
real::weight(newxy,4)
real newLon,newLat,weight_end  !本次插值的新点位置，以及返回的权重
real oldVal(oldx,oldy),newVal(newxy),newPoint(newxy,2)  !旧值，新值，新点经纬度
logical isfirst  !是否第一个时次，只有第一个时次需要计算四周格点分布以及权重，后面的套用即可

if(isfirst)then
!获取新点的四周的点以及每个点的权重
do i=1,newxy
  newLon=newPoint(i,1)

!print*,newLon,oldLon_start
  call getPosAndWeight(newLon,oldLon_start,nstart,nend,weight_end,oldInterpx,isFirst)
  !print*,newLon,nstart,nend,weight_end
  if(nend>oldx)then
    if(nstart==oldx)then
          pos(i,1)=nstart
          pos(i,2)=nstart
          weight(i,1)=1
          weight(i,2)=0
        else
      pos(i,1)=-1
          pos(i,2)=-1  !超出了原来格点的范围，标记为无效值
          weight(i,1)=-1
          weight(i,2)=-1
        endif
  else
    pos(i,1)=nstart
        pos(i,2)=nend
        weight(i,1)=1-weight_end
        weight(i,2)=weight_end
  endif

!print*, pos(i,1),pos(i,2)
!write(*,*) weight(i,1),weight(i,2)


  newLat=newPoint(i,2)
! write(*,*) newLat,oldLat_start
  call getPosAndWeight(newLat,oldLat_start,nstart,nend,weight_end,oldInterpy,isFirst)

  if(nend>oldy)then
    if(nstart==oldy)then
          pos(i,3)=nstart
          pos(i,4)=nstart
          weight(i,3)=1
          weight(i,4)=0
        else
      pos(i,3)=-1
          pos(i,4)=-1  !超出了原来格点的范围，标记为无效值
          weight(i,3)=-1
          weight(i,4)=-1
        endif
  else
    pos(i,3)=nstart
        pos(i,4)=nend
        weight(i,3)=1-weight_end
        weight(i,4)=weight_end
  endif

!print*, pos(i,3),pos(i,4)
!write(*,*) weight(i,3),weight(i,4)  
!pause   
enddo
endif



!开始插值
do i=1,newxy
  if(pos(i,1)<0 .or.pos(i,3)<0)then
    newval(i)=-999.0
  else
    !先做纬向插值
        if(pos(i,1)==undef .or. pos(i,2)==undef .or. pos(i,3)==undef .or. pos(i,4)==undef)then
          newVal(i)=-9999.0
        else
          if(oldVal(pos(i,1),pos(i,3))>0.and.oldVal(pos(i,2),pos(i,3))>0) tmpVal_down=weight(i,1)*oldVal(pos(i,1),pos(i,3))+weight(i,2)*oldVal(pos(i,2),pos(i,3))
          if(oldVal(pos(i,1),pos(i,3))>0.and.oldVal(pos(i,2),pos(i,3))<0) tmpVal_down=oldVal(pos(i,1),pos(i,3))
          if(oldVal(pos(i,1),pos(i,3))<0.and.oldVal(pos(i,2),pos(i,3))>0) tmpVal_down=oldVal(pos(i,2),pos(i,3))
          if(oldVal(pos(i,1),pos(i,3))<0.and.oldVal(pos(i,2),pos(i,3))<0) tmpVal_down=-999.0
          if(oldVal(pos(i,1),pos(i,4))>0.and.oldVal(pos(i,2),pos(i,4))>0) tmpVal_up=weight(i,1)*oldVal(pos(i,1),pos(i,4))+weight(i,2)*oldVal(pos(i,2),pos(i,4))
          if(oldVal(pos(i,1),pos(i,4))<0.and.oldVal(pos(i,2),pos(i,4))>0) tmpVal_up=oldVal(pos(i,2),pos(i,4))
      if(oldVal(pos(i,1),pos(i,4))>0.and.oldVal(pos(i,2),pos(i,4))<0) tmpVal_up=oldVal(pos(i,1),pos(i,4))
          if(oldVal(pos(i,1),pos(i,4))<0.and.oldVal(pos(i,2),pos(i,4))<0) tmpVal_up=-999.0
          if(tmpVal_down>0.and.tmpVal_up>0) then
                newVal(i)=weight(i,3)*tmpVal_down+weight(i,4)*tmpVal_up
          else if(tmpVal_down>0.and.tmpVal_up<0) then
            newVal(i)=tmpVal_down
      else if(tmpVal_down<0.and.tmpVal_up>0) then
                newVal(i)=tmpVal_up
          else
                newVal(i)=-999.0
          endif
        endif
  endif
enddo


open(25,file='E:\yh\weight.txt')
do i=1,newxy
write(25,*) weight(i,:)

enddo
close(25)

endsubroutine

!计算某个方向的两侧的点以及后侧点（上方点）的权重
subroutine  getPosAndWeight(newPos,oldstart,nstart,nend,weight,oldInterp,isfirst)
integer nstart,nend  !两侧点
real newPos,weight,oldInterp
real oldstart   
logical isfirst
!离的近的点权重大，但是占的距离比重小，所以权重=1-距离比重，经纬度从0开始，但是标号从1开始，所以加1
!print*,newPos,oldstart,oldInterp
nstart=int((newPos-oldstart)/oldInterp)+1
nend=nstart+1
weight=(newPos-(nstart-1)*oldInterp-oldstart)/oldInterp  !  |---|-------|  此表达式计算的是后一段的权重（前一段的距离比重），该段总长为interp
isfirst=.false.
endsubroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

站点信息数据：
50136 53.47  122.37
50353 51.72  126.65
50434 50.48  121.68
50442 50.40  124.12
50468 50.25  127.45
50514 49.57  117.43
50527 49.22  119.75
50557 49.17  125.23
50603 48.67  116.82