求助！！fortran读取753个站的逐日降水资料，并求月平均时遇到的问题

逸儿~~

这是老师给的“fortran读取全国753个站的逐日降水资料并求月平均降水”的程序（程序中的注释会说明资料的格式等），我稍作了下修改：

INTEGER,PARAMETER :: NS=756,NY=56+6,ND=31,NM=12,NY2=56
INTEGER year,mon,day
real    a(12),rainM(ns,ny,nm)
character*5 sta(NS)
real,allocatable ::  staID(:),lat(:),lon(:),hgt(:),pressure(:,:,:,:),temA(:,:,:,:),temH(:,:,:,:),temL(:,:,:,:), &
                  rh(:,:,:,:),rain(:,:,:,:),wind(:,:,:,:),sun(:,:,:,:)
allocate (staID(NS),lat(NS),lon(NS),hgt(NS),pressure(NS,NY,NM,ND),temA(NS,NY,NM,ND),temH(NS,NY,NM,ND),temL(NS,NY,NM,ND), &
      rh(NS,NY,NM,ND),rain(NS,NY,NM,ND),wind(NS,NY,NM,ND),sun(NS,NY,NM,ND))
!!!!!!日值说明，日平均气压、日平均气温、日最高气温、日最低气温、日平均相对湿度、日降水量、
!!!!!!          日平均风速、日日照时数
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  READ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
OPEN(11,FILE='F:\DATA\753-50-12\753-50-12\station.txt')
do i=1,NS
read(11,*)staID(i)
enddo
CLOSE(11)
  pressure=1.e36
  temA=1.e36
  temH=1.e36
  temL=1.e36
  rh  =1.e36
  rain=1.e36
  wind=1.e36
  sun =1.e36
!!  以下处理中将缺测、空白等值均定为1.e36，以方便处理
!===========1950-2005:
DO IS=1,NS
WRITE(sta(is),'(i5)')INT(staid(is))

OPEN(11,FILE='d:\DATA\753-50-12\753-50-12\DAY\SURF_CLI_CHN_MUL_DAY-'//sta(is)//'.txt',err=8000)
!!!!!!从老师的上述文件路径，可学到批处理许多文件的方法 !!!读取!!!
!!!!!!DAY文件夹里含753个站的txt资料，文件名都是“SURF_CLI_CHN_MUL_DAY-站号”的形式

!!!!!!!!区分读三个数据的程序的不同，注意数据的排放
200  read(11,*,err=100)(a(j),j=1,4),year,mon,day,(a(j),j=5,12)
        
  lat(is)=a(2)
  lon(is)=a(3)
  hgt(is)=a(4)
  pressure(is,year-1949,mon,day)=a(5)
  temA(is,year-1949,mon,day)=a(6)
  temH(is,year-1949,mon,day)=a(7)
  temL(is,year-1949,mon,day)=a(8)
  rh  (is,year-1949,mon,day)=a(9)
  rain(is,year-1949,mon,day)=a(10)
  wind(is,year-1949,mon,day)=a(11)
  sun (is,year-1949,mon,day)=a(12)
   goto 200
CLOSE(11)
100 continue
8000 continue
!!!!!  不是所有的台站都是从50年开始的，缺少的部分值均为1.e36，看下面的处理
!!!!!!200 与goto200 之间的程序表示 每读一天的数据，便将其赋值给日值数组，最终
!!!!!!将整个文件中的日值读取完
ENDDO
!!!!!!!整个do循环表示将文件夹D:\DATA\753-50-08\753-50-08\DAY中所有台站的资
!!!!!!!料都读取出来，赋值给日值数组pressure、temA、temH等、、、

!==========2006&2007:
DO IS=1,NS
WRITE(sta(is),'(i5)')INT(staid(is))
OPEN(11,FILE='d:\DATA\753-50-08\753-50-08\DAY\2006-07\SURF_CLI_CHN_MUL_DAY_0607-'//sta(is)//'.txt',err=400)
300  read(11,*,err=400)(a(j),j=1,4),year,mon,day,(a(j),j=5,12)
  pressure(is,year-1949,mon,day)=a(5)
  temA(is,year-1949,mon,day)=a(6)
  temH(is,year-1949,mon,day)=a(7)
  temL(is,year-1949,mon,day)=a(8)
  rh  (is,year-1949,mon,day)=a(9)
  rain(is,year-1949,mon,day)=a(10)
  wind(is,year-1949,mon,day)=a(11)
  sun (is,year-1949,mon,day)=a(12)
   goto 300
CLOSE(11)
400 continue
ENDDO
!!!!!!!此do循环程序和上面一个大致相同，只有一点，要注意50-05年有753个站，06-07年有728个站，所以open时，
!!!!!!!加err项，指定一个出错处理转移语句是必要的!!!
!!!!!!!50-07年，资料存放的格式是相同的，都是按照台站号、纬度、经度、位势高度、然后再是各日值，而08年
!!!!!!!的日值是按台站号、年、月、日、然后再是各日值摆放的，所以要注意它们的读取方式的不同

!==========2008:
OPEN(11,FILE='d:\DATA\753-50-08\753-50-08\08-11\day-08.txt')
     read(11,*)
!!!!!!read(11,*)表示从文件11中读出一个数据，但不写入任何变量
500  read(11,*,err=600)a(1:12)
  DO is=1,ns
  if(a(1).eq.staID(is))then
  year=a(2)
  mon =a(3)
  day =a(4)
  pressure(is,year-1949,mon,day)=a(5)
  temA(is,year-1949,mon,day)=a(6)
  temH(is,year-1949,mon,day)=a(7)
  temL(is,year-1949,mon,day)=a(8)
  rh  (is,year-1949,mon,day)=a(9)
  rain(is,year-1949,mon,day)=a(10)
  wind(is,year-1949,mon,day)=a(11)
  sun (is,year-1949,mon,day)=a(12)
  endif
        ENDDO
   goto 500
CLOSE(11)
!!!!!!!!因为前面日值数组都是按台站号、年月日排列的，所以对08年数据的处理也按同样方式，这就是do循
!!!!!环和if语句的作用
600 continue

!!!!!!!读完资料之后，日值都是756个站50-08年的，空白、缺测的部分均为1.e36
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  READ END  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!
!!!以下是对空白和缺测数据的处理

DO IS=1,NS
  lat(is)=INT(lat(is)/100)+MOD(lat(is),100.)/60.
  lon(is)=INT(lon(is)/100)+MOD(lon(is),100.)/60.
  if(hgt(is).ge.100000)hgt(is)=1.e36
!!!!!!当海拔高度为100000+xxxxx时，海拔高度为估测
!!  lat和lon是什么意思呢？这段程序是老师拷给的，可能和老师要做的内容有关，具体就不得而知了，我是这样想的
DO IY=1,NY
  do im=1,nm
  do id=1,nd
  if(abs(pressure(is,year-1949,mon,day)).eq.32744.or.abs(pressure(is,year-1949,mon,day)).eq.32766.or. &  
  pressure(is,year-1949,mon,day).eq.32700)then
  pressure(is,year-1949,mon,day)=1.e36
  endif
  if(abs(temA(is,year-1949,mon,day)).eq.32744.or.abs(temA(is,year-1949,mon,day)).eq.32766.or. &
  temA(is,year-1949,mon,day).eq.32700)then
  temA(is,year-1949,mon,day)=1.e36
  endif
  if(abs(temH(is,year-1949,mon,day)).eq.32744.or.abs(temH(is,year-1949,mon,day)).eq.32766.or. &
  temH(is,year-1949,mon,day).eq.32700)then
  temH(is,year-1949,mon,day)=1.e36
  endif
  if(abs(temL(is,year-1949,mon,day)).eq.32744.or.abs(temL(is,year-1949,mon,day)).eq.32766.or. &
  temL(is,year-1949,mon,day).eq.32700)then
  temL(is,year-1949,mon,day)=1.e36
!!!!!以上是对日平均气压、日平均气温、日最高气温和日最低气温在空白、缺测、微量或结冰情况下的处理均
!!!!!为1.e36，和初始值相同
  endif
  if(abs(rh(is,year-1949,mon,day)).eq.32744.or.abs(rh(is,year-1949,mon,day)).eq.32766)then
  rh(is,year-1949,mon,day)=1.e36
  elseif(rh(is,year-1949,mon,day).eq.32700)then
  rh(is,year-1949,mon,day)=0
  endif
!!!!!!以上是对日平均相对湿度的处理，空白、缺测情况下为1.e36，微量或结冰情况下为0
     if(abs(rain(is,year-1949,mon,day)).eq.32744.or.abs(rain(is,year-1949,mon,day)).eq.32766)then
  rain(is,year-1949,mon,day)=1.e36
  elseif(rain(is,year-1949,mon,day).eq.32700)then
  rain(is,year-1949,mon,day)=0
  elseif(rain(is,year-1949,mon,day).gt.30000.and.rain(is,year-1949,mon,day).lt.31000)then
  rain(is,year-1949,mon,day)=rain(is,year-1949,mon,day)-30000
  elseif(rain(is,year-1949,mon,day).gt.31000.and.rain(is,year-1949,mon,day).lt.32000)then
  rain(is,year-1949,mon,day)=rain(is,year-1949,mon,day)-31000
  elseif(rain(is,year-1949,mon,day).gt.32000.and.rain(is,year-1949,mon,day).lt.33000)then
  rain(is,year-1949,mon,day)=rain(is,year-1949,mon,day)-32000
  endif
!!!!!!以上是对日降水量的处理，空白、缺测情况下为1.e36，微量或结冰情况下为0，30xxx时，有雨雪，降
!!!!!!水量为xxx；31xxx时，有雪（雨夹雪、雪暴），降水量为xxx；32xxx时，有雾露霜，降水量为xxx
  if(abs(wind(is,year-1949,mon,day)).eq.32744.or.abs(wind(is,year-1949,mon,day)).eq.32766)then
  wind(is,year-1949,mon,day)=1.e36
  elseif(wind(is,year-1949,mon,day).gt.1000.and.wind(is,year-1949,mon,day).lt.2000)then
  wind(is,year-1949,mon,day)=wind(is,year-1949,mon,day)-1000
  endif
!!!!!!以上是对日平均风速的处理,空白、缺测情况下为1.e36，1000+xxx（即风速大于等于xxx时，加上1000）时，
!!!!!!风速为xxx
  if(abs(sun(is,year-1949,mon,day)).eq.32744.or.abs(sun(is,year-1949,mon,day)).eq.32766.or. &
  sun(is,year-1949,mon,day).eq.32700)then
  sun(is,year-1949,mon,day)=1.e36
  endif
!!!!!!以上是对日日照时数的处理,空白、缺测、微量或结冰情况下的处理均为1.e36
  
  enddo
  enddo
ENDDO
ENDDO
!===================== day=>month ================================================
!!!一下程序用于计算月平均降水量

DO IS=1,NS
DO IY=1,NY
do im=1,nm
   l=0
   rainM(is,iy,im)=0.
   do id=1,nd
   if(rain(is,iy,im,id).ne.1.e36)then
   l=l+1
!!!!!!! l在此处起计数作用，统计计算了多少天的降水
   rainM(is,iy,im)=rainM(is,iy,im)+rain(is,iy,im,id)
   endif
   enddo
   if((im.eq.1.or.im.eq.3.or.im.eq.5.or.im.eq.7.or.im.eq.8.or.im.eq.10.or.im.eq.12).and.l.eq.31)then
   rainM(is,iy,im)=rainM(is,iy,im)/l
!!  ???????这个地方感觉有问题啊，l只是统计降水不为1.e36的天数，而不是一个月有多少天，怎么判断条件里
!!!!!!!能有 .and.l.eq.31 呢??下面也是   是不是一个月里有缺测的，这个月的月平均降水就为缺测的呢？
   elseif((im.eq.4.or.im.eq.6.or.im.eq.9.or.im.eq.10.or.im.eq.11).and.l.eq.30)then
   rainM(is,iy,im)=rainM(is,iy,im)/l
   elseif(im.eq.2.and.l.ge.28)then
   rainM(is,iy,im)=rainM(is,iy,im)/l
   else
   rainM(is,iy,im)=1.e36
   endif
enddo
ENDDO
ENDDO
!=================================================================================
!do is=1,ns
!if(staID(is).eq.54511)ii=is
!!!!!! 54511是北京站
!enddo
!do is=1,ns
!if(  sqrt((lat(is)-lat(ii))**2+(lon(is)-lon(ii))**2).le.1.5)then
!print*,staid(is),is,lat(is),lon(is)
!!!!!! 找出离北京站较近的站点，并输出其台站号及经纬度
!endif
!enddo
!pause

OPEN(11,FILE='d:\data\756RainMon.dat')
DO IS=1,NS
DO IY=1,NY
  WRITE(11,*)staID(IS),hgt(IS),Lat(is),lon(is),(rainM(is,iy,im),im=1,12)
!!!!!!输出结果：中国756个站50-08年的月平均降水资料
ENDDO
ENDDO
CLOSE(11)
!!===============================================
deallocate (staID,lat,lon,hgt,pressure,temA,temH,temL,rh,rain,wind,sun)
END


程序编译之后是没有错的，但是运行后却出现一下结果：（运行后没有文件生成，亦不知错误在哪儿，还请大家指教）
下图中第二行，文件SURF_CLI_CHN_MUL_DAY-50136是DAY文件夹下要读取的05年的一个站的txt文件。。

做个霸气的木头

程序就不看了，实在长的吓人。一般出现这个问题，要么这个站号的文件没有，要么是这个文件的数据的排列不符合你的程序，导致文件无法正常读取，你需要带的是查对这个文件跟前面对的文件有什么区别，搞成一致的就好了。

river

首先要说read(11,*)表示的是空读一行，相当于跳过第一行，不是楼主说的那个。
再有l是计数的，但是你要看l在月循环里面被赋了初值0，所以表示的是每个月累加了多少天的降水。后面那个是判断大小月的吧。
最后那个错误说明SURF_CLI_CHN_MUL_DAY-50136资料长度比你设定的数组短

北五味

真不短啊，看着眼晕啊。

mofangbao

能不能把确定没错的自己先排除掉?

北五味

read里面不是一般都要加个end=endlabel以防读到文件末尾嘛。编程不熟，静待大神~~

逸儿~~

做个霸气的木头 发表于 2013-5-6 22:28

                               
登录/注册后可看大图

程序就不看了，实在长的吓人。一般出现这个问题，要么这个站号的文件没有，要么是这个文件的数据的排列不符 ...

不好意思啊，因为问过同学了都不知道错在哪儿了，不得已只能放在论坛上了，请论坛的高手指教，确实是长的吓人。。程序里面是有错误转移语句的，本来资料里面就是50-05年是753个站的，06-07年是728个站的，读取时也是按资料的排放读的。。出现错误的是读取的第一个文件。。

逸儿~~

river 发表于 2013-5-6 22:54

                               
登录/注册后可看大图

首先要说read(11,*)表示的是空读一行，相当于跳过第一行，不是楼主说的那个。
再有l是计数的，但是你要看l ...

嗯，谢谢啦，read那个问题虽然注释了，当时也不太确定的~~l的问题现在明白啦，如果l不等于那个月的天数的话，那个月的月平均就算缺测了~~再次感谢，嘿嘿~~

逸儿~~

北五味 发表于 2013-5-7 09:09

                               
登录/注册后可看大图

真不短啊，看着眼晕啊。

嘿嘿，不好意思哦，这个程序我也是断断续续看了好几天~~~

逸儿~~

mofangbao 发表于 2013-5-7 09:13

                               
登录/注册后可看大图

能不能把确定没错的自己先排除掉?

哎。。。我能说这个程序我已经检查了好几遍，弄了好几天了吗，我再一段一段的试一下~~谢谢学长~~