- 积分
- 3277
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-11-6
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
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
|
|
窥视卡
雷达卡
发表于 2018-4-19 15:16:32
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
{:5_217:}{:5_217:}{:5_217:}{:5_217:}{:5_217:}{:5_217:}