windows下EC细网格风场格点数据批量插值为站点数据求解

hushunqi215

在家园上看了帖子，想用fortran程序编写一个EC细网格风场资料的格点数据插值到站点数据的批量处理程序。主程序见附件

!   双线性插值，格点到站点，考虑权重  !
!      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

主程序编译没有错，运行后出现错误如下：run-time error F6200:READ(E:\yh\180522201.120)  
                                                             -formatted I/O not consistent with OPEN options
发此悬赏贴求解答。

hushunqi215

站点信息数据：
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

安达充

你读文件有问题

  open(29,file='E:\yh\18052220.120',form='formatted')
  DO i=1,5
    read(29,*)
  ENDDO
  DO j=1,oldy
    read(29,*)
    read(29,*) oldVal(:,j)
  ENDDO
  close(29)

hushunqi215

问题解决了，多谢