爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 3265|回复: 0

[求助] micaps 第二类数据转换为grads站点数据的问题

[复制链接]
发表于 2016-6-30 11:00:19 | 显示全部楼层 |阅读模式

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

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

x
各位大侠,
我最近写了一个fortran 程序(gfortran 编译器), 貌似一切都没有问题, 可是最后写成的grads站点数据,画图时, 400hp以上的数据都没有, 1000到400的数据貌似都是对的。 请各位大侠帮忙啊。现贴出我的程序和CTL文件。 现在的情况是: 程序正常运行,stnmap 也显示正常, 就是从300hpa以上就啥也显示不了了。

感觉就是300语句open或者write的不对呢

请大侠看看我open(300,然后的写法有什么错误么?  gfortran 写高空站点数据用什么格式好啊?  gfortran 不认binary

high_test.f90
!!提取高空的数据---一次处理一个月------------------
!!!!-------------------------------------------------------------------------------------------------
program main
implicit none
integer,parameter:: nn=800,lev=11,tt=2
!integer,parameter:: nn=800,lev=11, lev1=11,tt=2
real,parameter:: lo1=70,lo2=150,la1=15,la2=65
real,parameter::undef=9999
integer:: i,j,k,l,g,gg,m,it,j_s
integer:: nline(tt),num(tt)
integer:: nflag,nlev
integer:: ierr,irec
integer::kind1,kind2
real::tim,press(lev)                                                                                                   
real:: lon(nn,tt),lat(nn,tt),ls(nn,tt),hgt(nn,lev,tt),wd(nn,lev,tt),&
ttd(nn,lev,tt),fx(nn,lev,tt),fs(nn,lev,tt)
real:: lonh(nn,tt),lath(nn,tt),hgth(nn,lev,tt),wdh(nn,lev,tt),&
ttdh(nn,lev,tt),tdh(nn,lev,tt),fxh(nn,lev,tt),fsh(nn,lev,tt)
real:: est(nn,lev,tt), estd(nn,lev,tt), qvs(nn,lev,tt),qes(nn,lev,tt)
real :: rh(nn,lev,tt)   
integer:: zj(nn,tt),sta(nn,tt),stah(nn,tt),st(nn)
character(len=2)::yy,hour
character(len=4)::yr
character(len=2)::mon
character(len=100)::dirname,data_in,data_out
character(len=200)::filename, openname
character:: a(4)*1,b(10)*1
character(len=4):: hp(lev)
character(len=2):: sh(2)
character*8 stid(nn)
logical alive,live
character(len=10):: temp(5)
data a/'0','1','2','3'/
data b/'0','1','2','3','4','5','6','7','8','9'/
data hp/'1000','925','850','700','500','400','300','250','200','150', &
        '100'/
data sh/'08','20'/
REAL    , PARAMETER ::  SVP1=0.6112
REAL    , PARAMETER ::  SVP2=17.67
REAL    , PARAMETER ::  SVP3=29.65
REAL    , PARAMETER ::  SVPT0=273.15
REAL :: ALIQ, BLIQ, CCLIQ, DLIQ
integer:: int_yr,start_date, end_date
integer:: dates(12)=(/ 31,28,31,30,31,30,31,31,30,31,30,31/)
integer:: reck
!-----------------------------屏幕读入处理时间-------------------------------
write(*,*) "enter year month:******(201401)"
read(*,"(a4,a2)")yr,mon
yy=yr(3:4)
read(yr, "(i4)") int_yr
write(*,*) "enter start_date:(01)"
read(*,"(i2)")start_date
write(*,*) "enter end_date:(31)"
read(*,"(i2)")end_date
write(*,*) "enter finish"
write(*,*) "yr,mon,start_date,end_date", yr,mon,start_date,end_date
!--------------------------------------判断是否闰年----------------------------------------
if((mod(int_yr,4).eq.0.and.mod(int_yr,100).ne.0) .or. mod(int_yr,400).eq.0) then
   dates(2)=29
else
   dates(2)=28
endif
!------------------------------------prepare for rh calculating -----------------------------
  ALIQ = SVP1*1000.
     BLIQ = SVP2
     CCLIQ = SVP2*SVPT0
     DLIQ = SVP3
!!!---------------------输入、输出文件夹---------------------------
data_in="D:\obs_pro\data\micaps\"
data_out="D:\obs_pro\output\met_data\high\high_stn\"
do l=1,lev
read(hp(l),*) press(l)
end do
write(*,*) press(:)
!!!!----
do j=start_date, end_dateme,'(3a,i2.2,2a)') trim(data_out),trim(yr),trim(mon),j,sh(k),'.grd'
    write(*,*) 'output name', filename
!   OPEN(300,FILE=trim(filename),FORM='unformatted', recordtype=' stream')
OPEN(300,FILE=trim(filename),FORM='unformatted',access='sequential')
do l=1,lev            
   write(openname,'(8a,i2.2,2a)') trim(data_in),trim(yr),trim(mon),'\high\plot\',trim(adjustl(hp(l))), &
    '\',trim(yy),trim(mon),j,sh(k),".000"
   print*,'read file', openname
   open(100,file=trim(openname),status='old')
      read(100,*)
      read(100,*) (temp(it),it=1,5),nline(k)
   print*, 'nline(k)', nline(k)
    do g=1,nline(k)
      read(100,*)sta(g,k),lon(g,k),lat(g,k),ls(g,k),zj(g,k),hgt(g,l,k),wd(g,l,k),&
               ttd(g,l,k),fx(g,l,k),fs(g,l,k)      
end do !g
    close(100)
  end  do
  

do m=1,nline(k)
   write(stid(m),'(i8)') stah(m,k)  ! convert integer to character
enddo !m
  
print*, 'before write 300 lev='
  tim=0.0
  nlev=lev
  nflag=0
  irec=1
  do m=1,nline(k)
      
    write(300)  stid(m), lat(m,k),lon(m,k),tim,nlev,nflag
      
    !    write(300)  stid(m), lat(m,k),lon(m,k),tim,5,nflag
do l=1,lev
!    do l=1,5
    print*, 'before write 300, st,lev,press', m,l
write(300)  press(l),hgt(m,l,k),wd(m,l,k),ttd(m,l,k),fx(m,l,k),fs(m,l,k)

print*, press(l),hgt(m,l,k),wd(m,l,k),ttd(m,l,k),fx(m,l,k),fs(m,l,k)  
print*, 'press,hgt,wd(m,l,k),ttd(m,l,k),fx(m,l,k),fs(m,l,k)'

enddo !l1
  enddo !m
  
  nlev=0
       write(300)  stid(m-1), lat(m-1,k),lon(m-1,k),tim,nlev,nflag
   
close(300)
end do  
end do  
end program

下面是CTL文件high_test.ctl
dset D:\obs_pro\output\met_data\high\high_stn\2015110120.grd
dtype station
options sequential
stnmap  D:\obs_pro\output\met_data\high\high_stn\2015110120.map
undef  9999
title high
TDEF 1  linear  12Z01NOV2015  12hr
vars 5
hgt  11 99 geopotential height
temp 11 99 temperature
ttd 11 99 T-Td
windd  11 99 wind direction
winds  11 99 wind speed
endvars


密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

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

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

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