[已解决]753站点提取数据画降水图多时次问题

echo要读书

我想在573站中提取出1991-2005年的5-8月逐日降水资料。问题：映射文件出错－－－》
程序如下
program main
integer*4,parameter::ns=753,ny=56,nd=31,nm=12
integer*4::year,mon,day,nlev,flag,is
real*4::a(12),tim
integer*4::mdays(12)
character*8::sta(ns)
!character*8,allocatable::staid(:)
real,allocatable::lat(:),lon(:),hgt(:),pressure(:,:,:,:),temA(:,:,:,:),temH(:,:,:,:),temL(:,:,:,:), &
                  rh(:,:,:,:),rain(:,:,:,:),wind(:,:,:,:),sun(:,:,:,:)
allocate (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 station names!!!!!!!!!!!!!!!!!!!!!

open(11,file='z:\prog\grd\fro\7535012\station.txt')
do i=1,ns
read(11,*) sta(i)
write(*,*) sta(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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!read data!!!!!!!
do is=1,ns

open (11,file='Z:/prog/grd/fro/7535012/7535012/SURF_CLI_CHN_MUL_DAY-'//trim(adjustl(sta(is)))//'.txt',err=8000)

200 read(11,*,end=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

enddo


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!wrong data!!!!!!!!!!!!!!!!!!!!!

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

  do iy=1,ny
     do im=1,nm
            do id=1,nd
                   if(abs(rain(is,iy,im,id)).eq.32744.or.abs(rain(is,iy,im,id)).eq.32766)then
             rain(is,iy,im,id)=1.e36
            elseif(rain(is,iy,im,id).eq.32700)then
             rain(is,iy,im,id)=0
                        elseif(rain(is,iy,im,id).gt.30000.and.rain(is,iy,im,id).lt.31000)then
                         rain(is,iy,im,id)=rain(is,iy,im,id)-30000
                        elseif(rain(is,iy,im,id).gt.31000.and.rain(is,iy,im,id).lt.32000)then
                         rain(is,iy,im,id)=rain(is,iy,im,id)-31000
                        elseif(rain(is,iy,im,id).gt.32000.and.rain(is,iy,im,id).lt.33000)then
                         rain(is,iy,im,id)=rain(is,iy,im,id)-32000
                        endif
                    enddo
                enddo
        enddo                                                                                                                                                                                 
enddo

!!!!!!!!!!!!!!!!!!!!!!!write!!!!!!!!!!!!!!!!!!!!!!!!!

mdays(5) = 31
mdays(6) = 30
mdays(7) = 31
mdays(8) = 31

OPEN(11,FILE='Z:\prog\grd\fro\7535012\756precip.grd')
DO IY=42,NY

     Do IM=5,8
      write(*,*) im
      do id=1,mdays(im)
      tim=0.0
      nlev=1
      flag=1
      do is=1,ns
      WRITE(11,*)sta(IS),Lat(is),lon(is),tim,nlev,flag,rain(is,iy,im,id)
          WRITE(*,*)sta(IS),iy,im,id,Lat(is),lon(is),tim,nlev,flag,rain(is,iy,im,id)
      enddo
      nlev=0
      WRITE(11,*)sta,Lat,lon,tim,nlev,flag
      enddo
        enddo
enddo
close(11)

!!!!!!!!==================================
deallocate (lat,lon,hgt,pressure,temA,temH,temL,rh,rain,wind,sun)
END

river

grads要求站点资料写入方式就是每个时次必须有结束符，无论单时次还是多时次都需要。还有站点stid 必须是 character*8的，这个是特定的不许改。你不是看了清风的帖子吗，怎么把最主要的东西都给改了呢······

echo要读书

river 发表于 2014-7-3 21:22
grads要求站点资料写入方式就是每个时次必须有结束符，无论单时次还是多时次都需要。还有站点stid 必须是 c ...

WRITE(sta(is),'(i8)')INT(staid(is))
＝
没有改啊？

river

echo要读书 发表于 2014-7-3 22:17
WRITE(sta(is),'(i8)')INT(staid(is))
＝
没有改啊？

那就忽略后一句吧······

lqouc

echo要读书 发表于 2014-7-3 22:17
WRITE(sta(is),'(i8)')INT(staid(is))
＝
没有改啊？

是么？那你输出的是staid还是sta？
不管怎么样，人家帮你解决问题，不要急着反驳，试一下他的建议，如果正确就表示感谢。

echo要读书

lqouc 发表于 2014-7-3 22:35
是么？那你输出的是staid还是sta？
不管怎么样，人家帮你解决问题，不要急着反驳，试一下他的建议，如果 ...

额我没反驳啊。。这真的是问句。。很久没用fortran了，都是临时对着教程学的，所以对那个语句的正确性不太自信，数据在办公室的电脑里面现在没法试，才想听river进一步的解释一下的，真心毫无攻击性{:sweat:}{:sweat:}……我理解的这个语句是把staid（is）以8位int的形式写进sta(is)，这样不能将动态数组也变成8位么？你们的意思是最开始就要先把staid（：）定义成character*8, allocatable staid(:)这样?..
可能我的问题太小白了＝ ＝,还请你们指点一下... ...

echo要读书

river 发表于 2014-7-3 22:25
那就忽略后一句吧······

见楼上。。。otz如果我上层回复看起来语气有问题我道歉哈。。。。。。。= =!

river

echo要读书 发表于 2014-7-3 23:28
见楼上。。。otz如果我上层回复看起来语气有问题我道歉哈。。。。。。。= =!

这个和fortran是没有关系的，这种格式是grads规定的，他只能识别这样的站点数据。除此站点格式外，grads都是无法识别的。So，staid必须在变量生声明的时候就要定义成character*8的，或者变量声明的时候定义一个别的变量是character*8，最后写入的时候让这个变量等于后面处理的staid。

最后还想说你回复了几次，每次都纠结站点这个character*8，似乎是忽略了我最前面的一句话······但是那也是grads能识别的站点资料格式的一个要求，那一点达不到，生成stnmap的时候还是要出错······

echo要读书

river 发表于 2014-7-4 06:08
这个和fortran是没有关系的，这种格式是grads规定的，他只能识别这样的站点数据。除此站点格式外，grads ...

根据您的建议又修改了好几次。。。还是不能生成MAP文件。。泪奔TAT

echo要读书

river 发表于 2014-7-4 06:08
这个和fortran是没有关系的，这种格式是grads规定的，他只能识别这样的站点数据。除此站点格式外，grads ...

新修改的程序如下，省去缺省值处理的部分。目的是把1991-2005年5-8月的数据提取出来改成二进制
可以生成二进制文件，但是制作MAP文件时报错依然如主楼。

program main
integer*4,parameter::ns=753,ny=56,nd=31,nm=12
integer*4::year,mon,day,nlev,flag,is
real*4::a(12),tim
integer*4::mdays(12)
character*8::sta(ns)
!character*8,allocatable::staid(:)
real,allocatable::lat(:),lon(:),hgt(:),pressure(:,:,:,:),temA(:,:,:,:),temH(:,:,:,:),temL(:,:,:,:), &
                  rh(:,:,:,:),rain(:,:,:,:),wind(:,:,:,:),sun(:,:,:,:)
allocate (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 station names!!!!!!!!!!!!!!!!!!!!!

open(11,file='z:\prog\grd\fro\7535012\station.txt')
do i=1,ns
read(11,*) sta(i)
write(*,*) sta(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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!read data!!!!!!!
do is=1,ns

open (11,file='Z:/prog/grd/fro/7535012/7535012/SURF_CLI_CHN_MUL_DAY-'//trim(adjustl(sta(is)))//'.txt',err=8000)

200 read(11,*,end=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

enddo


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!wrong data!!!!!!!!!!!!!!!!!!!!!

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

  do iy=1,ny
     do im=1,nm
            do id=1,nd
                   if(abs(rain(is,iy,im,id)).eq.32744.or.abs(rain(is,iy,im,id)).eq.32766)then
             rain(is,iy,im,id)=1.e36
            elseif(rain(is,iy,im,id).eq.32700)then
             rain(is,iy,im,id)=0
                        elseif(rain(is,iy,im,id).gt.30000.and.rain(is,iy,im,id).lt.31000)then
                         rain(is,iy,im,id)=rain(is,iy,im,id)-30000
                        elseif(rain(is,iy,im,id).gt.31000.and.rain(is,iy,im,id).lt.32000)then
                         rain(is,iy,im,id)=rain(is,iy,im,id)-31000
                        elseif(rain(is,iy,im,id).gt.32000.and.rain(is,iy,im,id).lt.33000)then
                         rain(is,iy,im,id)=rain(is,iy,im,id)-32000
                        endif
                    enddo
                enddo
        enddo                                                                                                                                                                                 
enddo

!!!!!!!!!!!!!!!!!!!!!!!write!!!!!!!!!!!!!!!!!!!!!!!!!

mdays(5) = 31
mdays(6) = 30
mdays(7) = 31
mdays(8) = 31

OPEN(11,FILE='Z:\prog\grd\fro\7535012\756precip.grd')
DO IY=42,NY

     Do IM=5,8
      write(*,*) im
      do id=1,mdays(im)
      tim=0.0
      nlev=1
      flag=1
      do is=1,ns
      WRITE(11,*)sta(IS),Lat(is),lon(is),tim,nlev,flag,rain(is,iy,im,id)
          WRITE(*,*)sta(IS),iy,im,id,Lat(is),lon(is),tim,nlev,flag,rain(is,iy,im,id)
      enddo
      nlev=0
      WRITE(11,*)sta,Lat,lon,tim,nlev,flag
      enddo
        enddo
enddo
close(11)

!!!!!!!!==================================
deallocate (lat,lon,hgt,pressure,temA,temH,temL,rh,rain,wind,sun)
END