请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 9133|回复: 6

[源代码] 给大家分享一个格点插值到站点然后TS评分的程序

[复制链接]

新浪微博达人勋

发表于 2018-4-19 15:16:32 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
!****************************************************************************
!
!  PROGRAM: T639格点资料插值成站点资料
!
!****************************************************************************

    module pra_mod
           implicit none
           integer,parameter :: Nlon=90                               ! T639 E0-180  2
           integer,parameter :: Nlat=45                               ! T639 N0-90   2
           integer,parameter :: Nsta_max=5000                         ! 最多5000站                     
       real,parameter :: yuezhi=0.1                               ! 有无降水的阈值  mm
           character(len = *), parameter :: dir_T639='C:\Users\Tianle\Desktop\插值检验\t639-05\'        ! input original_data
           character(len = *), parameter :: file1_t639_name='16112102.027'       ! input original_data 体现那一天
           character(len = *), parameter :: file2_t639_name='16112202.027'       ! input original_data
       character(len = *), parameter :: dir_station='C:\Users\Tianle\Desktop\插值检验\'             ! input original_data
       character(len = *), parameter :: filename_station='STATIONS.DAT'      ! input original_data
       character(len = *), parameter :: dir_obs='C:\Users\Tianle\Desktop\插值检验\sur-05-24\'       ! input original_data
              character(len = *), parameter :: file1_obs_name='16112205.000'        ! input original_data 文件日期
           character(len = *), parameter :: file2_obs_name='16112305.000'        ! input original_data

        endmodule pra_mod

        module data_mod
           use pra_mod
           implicit none
           real lon_t639(0:Nlon)
           real lat_t639(0:Nlat)
           real Rain_t639(0:Nlon,0:Nlat)
           integer Nsta_t639,Nsta_obs,Nsta_common
           integer sta_number_t639(Nsta_max)
           real sta_lat_t639(Nsta_max),sta_lon_t639(Nsta_max),Rain_t639_sta(Nsta_max)
       character(len=10) sta_name_t639(Nsta_max)
           real Rain_obs(Nsta_max)
           integer sta_number_obs(Nsta_max)
           integer sta_number_common(Nsta_max)
           real Rain_t639_common(Nsta_max),Rain_obs_common(Nsta_max)




    end module data_mod

         
        program data_chazhi   ! main program
           use pra_mod
           use data_mod
       implicit none
       character(len=50)  cha
           character(len=10)  cha10
           integer station,stat,templat,templon,temp1,temp2
           real lon,lat,gh,rain,ts
           real  rf,lf,uf,df,lonP,latP
           integer rn,ln,un,dn,i,j

           do i=0,Nlon
              lon_t639(i)=i*2
                  !write(*,*) "lon",i,lon_t639(i)
           enddo

           do j=0,Nlat
              lat_t639(j)=j*2
                  !write(*,*) "lat",j,lat_t639(j)
           enddo

       write(*,*)'step1'
       Rain_t639(:,:)=0.
       open(22,file=dir_T639//file1_t639_name,form="formatted")
          
       read(22,*) cha
           !write(*,*) cha
           read(22,*) cha
           !write(*,*) cha
           do while(1)
          read(22,'(i6,f8.3,f8.3,f8.1,x,f8.1)',iostat=stat)   station,lon,lat,gh,rain
                  if (stat.lt.0)  exit  ! stat.lt.0  文件结束
          i=lon/2
                  j=lat/2
                  !write(*,*) i,j,rain
          Rain_t639(i,j)=rain
           enddo
           close(22)
           i=90
           j=45
       write(*,'(f8.1,f8.1,4x,f8.1)')   i*2.0,j*2.0,Rain_t639(i,j)

       Nsta_t639=0
       open(22,file=dir_station//filename_station,form="formatted")
           do while(1)
          read(22,*,iostat=stat)   station,templat,templon,gh,temp1,temp2,cha10
                  if (stat.lt.0)  exit  ! stat.lt.0  文件结束
                  Nsta_t639=Nsta_t639+1
          !write(*,*) Nsta_t639,station,templat,templon,cha10
          sta_number_t639(Nsta_t639)=station
          sta_lat_t639(Nsta_t639)=templat/100.
          sta_lon_t639(Nsta_t639)=templon/100.
                  sta_name_t639(Nsta_t639)=cha10         
           enddo
           close(22)
           write(*,'(a20,i6)')   "stations number",Nsta_t639

       ! test integrate_point
       !lonP=115.7
           !latP=36.2
       !call integrate_point(lon_t639,Nlon,lat_t639,Nlat,lonP,latP, &
           !                         rf,lf,uf,df,rn,ln,un,dn)


           open(22,file=dir_station//"Rain_t639_station.txt",form="formatted")
           do i=1,Nsta_t639
          lonP=sta_lon_t639(i)
              latP=sta_lat_t639(i)
          call integrate_point(lon_t639,Nlon,lat_t639,Nlat,lonP,latP, &
                                    rf,lf,uf,df,rn,ln,un,dn)
                  Rain_t639_sta(i)= Rain_t639(rn,un)*rf*uf &
                              +Rain_t639(rn,dn)*rf*df &
                                          +Rain_t639(ln,dn)*lf*df &
                      +Rain_t639(ln,un)*lf*uf
          if (((rn+ln).eq.0).or.((un+dn).eq.0)) then
                     write(*,'(a10,i5,2x,2f5.1)')  "bianjie:",i,sta_lon_t639(i),sta_lat_t639(i)   ! 超出边界
             sta_number_t639(i)=0
                  endif
                  write(22,'(i5,2x,2f5.1,2x,4f6.1,x,4f7.2,2x,i6,f7.2)')  i,sta_lon_t639(i),sta_lat_t639(i), &
                      lon_t639(ln),lon_t639(rn),lat_t639(un),lat_t639(dn), &
                      Rain_t639(rn,un),Rain_t639(rn,dn),Rain_t639(ln,un),Rain_t639(ln,dn),sta_number_t639(i),Rain_t639_sta(i)

           enddo
           close(22)



       Rain_obs(:)=0
           open(22,file=dir_obs//file1_obs_name,form="formatted")!jianyandiertianjiangshui
           do i=1,14
          read(22,*) cha
              !write(*,*) cha
       enddo
           Nsta_obs=0
           do while(1)
          read(22,'(i6,f8.3,f8.3,f6.0,f6.1)',iostat=stat)   station,lon,lat,gh,rain
                  if (stat.lt.0)  exit  ! stat.lt.0  文件结束
                  Nsta_obs=Nsta_obs+1
          sta_number_obs(Nsta_obs)=station
          Rain_obs(Nsta_obs)=rain
                  !write(*,*)  Nsta_obs,sta_number_obs(Nsta_obs), Rain_obs(Nsta_obs)
           enddo
           close(22)
       write(*,'(a20,i6)')   "stations_obs number",Nsta_obs


       open(22,file=dir_station//"Rain_t639_obs_common.txt",form="formatted")
       Nsta_common=0
       do i=1,Nsta_t639
          stat=0
                  do j=1,Nsta_obs
             if (sta_number_t639(i).eq.sta_number_obs(j)) then
                            stat=1
                                exit
             endif
                  enddo
                    Nsta_common = i
                  if (stat.eq.0)  then
!                     sta_number_t639(i)=0
                 write(22,'(i5,2i6,2x,2f7.1)')i,sta_number_t639(i),sta_number_t639(i),Rain_t639_sta(i),0
            Rain_t639_common(Nsta_common)=Rain_t639_sta(i)
                         Rain_obs_common(Nsta_common)=0
                  else

             Rain_t639_common(Nsta_common)=Rain_t639_sta(i)
                         Rain_obs_common(Nsta_common)=Rain_obs(j)
                         write(22,'(i5,2i6,2x,2f7.1)') i,sta_number_t639(i),sta_number_obs(j),Rain_t639_sta(i),Rain_obs(j)
                  endif
           enddo

           write(*,'(a20,i6)')   "Nsta_common number",Nsta_common



       call QETS(Nsta_common,Rain_t639_common(1:Nsta_common),Rain_obs_common(1:Nsta_common),ts)
       open(33,file='C:\Users\Tianle\Desktop\插值检验\ts\ts1.txt')

       write(33,*) "ts:",ts



          

        endprogram data_chazhi



        subroutine QETS(N,forc,obse,ts)  
          use pra_mod,only:yuezhi
      implicit none
          integer :: N                ! 用于评分的站点数目
      real :: forc(N),obse(N)     ! 预报和观测的降水序列
      real :: ts                  ! ts得分
          integer :: A,B,C            ! 命中,空报,漏报
          integer :: i
          A=0
          B=0
          C=0
          write(*,*)'yuezhi',yuezhi
      do i=1,N
         if(forc(i)>=yuezhi.and.obse(i)>=yuezhi) A=A+1
         if(forc(i)>=yuezhi.and.obse(i)<yuezhi)  B=B+1
         if(forc(i)<yuezhi.and.obse(i)>=yuezhi)  C=C+1
      enddo
      ts=A*1.0/(A+B+C)
      write(*,*) A,B,C
    end subroutine QETS


        subroutine integrate_point(lonMap,Ncols,latMap,Nrows,lonP,latP, &
                                    rf,lf,uf,df,rn,ln,un,dn)
           implicit none
           integer,intent(in) :: Ncols,Nrows
           real,intent(in) :: lonMap(0:Ncols),latMap(0:Nrows)
           real,intent(in) :: lonP,latP
           real,intent(out) :: rf,lf,uf,df
           integer,intent(out) :: rn,ln,un,dn
           integer i,j
           if ((lonP.le.lonMap(0)).or.(lonP.ge.lonMap(Ncols))) then   ! 超出边界 rn=ln=0
          rn=0
                  ln=0
                  rf=0.
                  lf=0.
                  write(*,*) " (lonP.le.lonMap(0)).or.(lonP.ge.lonMap(Ncols))",lonP
                  return
           endif
           if ((latP.le.latMap(0)).or.(latP.ge.latMap(Nrows))) then   ! 超出边界  un=dn=0
          un=0
                  dn=0
                  uf=0.
                  df=0.
                  write(*,*) "(latP.le.latMap(0)).or.(latP.ge.latMap(Nrows))",latP
                  return
           endif
       rn=1
       ln=1
       rf=1.
       lf=0.
       do i=1,Ncols
          if (lonP.le.lonMap(i)) exit
       enddo
       rn=i
       ln=i-1
       rf=(lonP-lonMap(i-1))/(lonMap(i)-lonMap(i-1))
       lf=(lonMap(i)-lonP)/(lonMap(i)-lonMap(i-1))            
           !write(*,'(2i5,6f8.3)') ln,rn,lf,rf,lonP,lonMap(ln),lonMap(rn),lf+rf
       un=1
       dn=1
       uf=1.
       df=0.
       do j=1,Nrows
          if (latP.le.latMap(j)) exit
       enddo
       un=j
       dn=j-1
       uf=(latP-latMap(j-1))/(latMap(j)-latMap(j-1))
       df=(latMap(j)-latP)/(latMap(j)-latMap(j-1))  
       !write(*,'(2i5,6f8.3)') dn,un,df,uf,latP,latMap(dn),latMap(un),df+uf                
    end subroutine integrate_point








密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-7-16 16:28:15 | 显示全部楼层
这个里面用的什么插值方法呀?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-7-25 16:36:15 | 显示全部楼层
谢谢分享!!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2018-8-9 10:00:59 | 显示全部楼层
请问用的什么插值方法
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-5-21 14:54:53 | 显示全部楼层
谢谢分享,学习了!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-3-31 14:05:14 | 显示全部楼层
{:5_217:}{:5_217:}{:5_217:}{:5_217:}{:5_217:}{:5_217:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-12-5 10:23:55 | 显示全部楼层

谢谢楼主无私分享!!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表