【推荐】站点降水数据批量处理程序（最近使用的）

kongfeng0824

以下是自己最近处理站点数据的程序，包括了闰年和平年，大月和小月的判别。分享给大家
批量处理降水数据：


program main
!integer,parameter:: years=2009,yearo=2009,months=07,montho=07,days=01,dayo=01,ipar=05,istations=2000,nx=1200,ny=600
include 'parrf.cb'
integer id,iidt(ipar),iid(ipar),monthz(12),irfav(istations),&
        iyear,imonth,month(12),monthb,monthf,dayb,dayf,iday,rn,nstation(31),position(31),ipart,idata,istation,&
                nlev,nflag,anlev,anflag,avnlev(ipar,istations),avnflag(ipar,istations),ia

real*4 grid(nx,ny),lon(istations,31),lat(istations,31),rf(istations,31),rfr(istations),hgt,tim,&
       avlat(ipar,istations),avlon(ipar,istations),avtime(ipar,istations),&
           avrf(ipar,istations),alat,alon,atime,rfrav(istations)
character ctime(31)*6,stid*8,astid*8,avstid(ipar,istations)*8,ro*10,name(istations,31)*5,cpar*2,&
          path_pic*9,cpath_pic*100,pathin*7,pathout*8,cpathin*100,cpathout*100,pathriver*8,cpathriver*100
logical alive
!----------
!head record format need vars
!character stid*8
!real lat
!real lon
!real tim
!integer nlev
!integer nflag
!------------------------
!read parrf.txt
open(1,file='.\parrf.txt')
read(1,*)pathin,cpathin
read(1,*)pathout,cpathout
read(1,*)pathriver,cpathriver
read(1,*)path_pic,cpath_pic
close(1)
!------------------------
!create gridfile for using
!china.grd
do i=1,nx
do j=1,ny
  grid(i,j)=1.0
enddo !j
enddo !i
open(1,file=cpathriver(1:len_trim(cpathriver))//'china.grd',Access='Direct',Form = 'Unformatted',recl=nx*ny)
write(1,rec=1)((grid(i,j),i=1,nx),j=1,ny)
close(1)
!china.ctl
open(1,file=cpathriver(1:len_trim(cpathriver))//'china.ctl')
write(1,'(A)')'dset  '//cpathriver(1:len_trim(cpathriver))//'china.grd'
write(1,'(A)')'title Sample GRIB Data'
write(1,'(A)')'undef -999.0'
write(1,'(A)')'xdef 1200 linear 60.0 0.1'
write(1,'(A)')'ydef 600 linear 0.0 0.1'
write(1,'(A)')'zdef 1 levels 1000'
write(1,'(A)')'tdef 1 linear 18z19Aug2009 1dy'
write(1,'(A)')'vars  1'
write(1,'(A)')'g   0  99  grid data prepared for oacres function'
write(1,'(A)')'endvars'
close(1)

!--------------------------
write(cpar,'(2I2.2)')ipar
do iyear=years,yearo
rn=iyear-iyear/4*4
if(rn.eq.0) month=[31,29,31,30,31,30,31,31,30,31,30,31]
if(rn.ne.0) month=[31,28,31,30,31,30,31,31,30,31,30,31]

if(iyear.eq.years) then
  monthb=months
else
  monthb=1
endif !monthA
if(iyear.eq.yearo) then
  monthf=montho
else
  monthf=12
endif !monthB

do imonth=monthb,monthf
  if((iyear.eq.years).and.(imonth.eq.months)) then
   dayb=days
  else
   dayb=01
  endif !dayA
  if((iyear.eq.yearo).and.(imonth.eq.montho)) then
   dayf=dayo
  else
   dayf=month(imonth)
  endif !dayB

  do iday=dayb,dayf
!-------------------------------------
!write name   
   write(ctime(31),'(3I2.2)')iyear-iyear/100*100,imonth,iday
   do it=1,30
    if(imonth.eq.1) then !A1
     if(it.lt.iday) then !A1B1
      write(ctime(it),'(3I2.2)')iyear-iyear/100*100,imonth,iday-it
     else !A1B2
      write(ctime(it),'(3I2.2)')(iyear-1)-(iyear-1)/100*100,imonth-1+12,iday-it+month(imonth-1)
     endif !A1B3
    else !A2
         if(it.lt.iday) then !A2B1
      write(ctime(it),'(3I2.2)')iyear-iyear/100*100,imonth,iday-it
     else !A2B3
          if(iday-it+month(imonth-1).le.0) write(ctime(it),'(3I2.2)')iyear-iyear/100*100,imonth-2,iday-it+month(imonth-1)+month(imonth-2)
          if(iday-it+month(imonth-1).gt.0) write(ctime(it),'(3I2.2)')iyear-iyear/100*100,imonth-1,iday-it+month(imonth-1)
     endif !A2B3
        endif !A3
   enddo !it
   print*,ctime(31)

   do i=1,31
    inquire(file=cpathin(1:len_trim(cpathin))//ctime(i)//'08.000',exist=alive)
    if(alive) position(i)=1
    if(alive.eq..False.) position(i)=0
   enddo
   !print*,position

   !-------------------
   !check files
   if(position(31).eq.0) print*,'data ',ctime(31),' do not exsit'
   if(position(31).eq.1) then !data exsit do the read
    i=1
    open(1,file=cpathin(1:len_trim(cpathin))//ctime(31)//'08.000')
    do k=1,14
     read(1,*)ro
!   print*,ro
    enddo !k
    do while(.true.)
     read(1,*,iostat=io) name(i,31),lon(i,31),lat(i,31),hgt,rf(i,31)
     if(io.lt.0) goto 100
         i=i+1
    enddo !while
    100 close(1)
    nstation(31)=i-1
       
    !---------------
    !check data before
    do ipart=1,ipar
     if(position(ipart).eq.0) print*,'data ',ctime(ipart),' do not exsit'
         if(position(ipart).eq.1) then
      i=1
      open(1,file=cpathin(1:len_trim(cpathin))//ctime(ipart)//'08.000')
      do k=1,14
       read(1,*)ro
      enddo !k
      do while(.true.)
       read(1,*,iostat=io) name(i,ipart),lon(i,ipart),lat(i,ipart),hgt,rf(i,ipart)
       if(io.lt.0) goto 200
           i=i+1
      enddo !while
      200 close(1)
     nstation(ipart)=i-1
        ! print*,nstation
         endif !position
    enddo !ipart
    !-----------------
    !do the math
    do istation=1,nstation(31)
     i=0
         do ipart=1,ipar
      do ist=1,nstation(ipart)
           if(name(istation,31).eq.name(ist,ipart)) then
            rfr(istation)=rfr(istation)+rf(ist,ipart)
            i=i+1
           endif !name
          enddo !ist
         enddo !ipart
         if(i.eq.0) then
          rfr(istation)=-999.0
!         else
!     rfr(istation)=rfr(istation)/i
         endif !i
    enddo !istation
!        print*,rfr(1:nstation(31))


  !=============================================================
  !------------------
        !read history av_data
        !confirm the day
        monthz=[31,29,31,30,31,30,31,31,30,31,30,31]
        id=0
        iidt=0
    if(imonth.eq.1) then
         id=iday
        else
         do i=1,imonth-1
          id=id+monthz(i)
         enddo !i
         id=id+iday
        endif !imonth
        do i=1,ipar
         if(id-i.le.0) then
          iidt(i)=id-i+366
         else
          iidt(i)=id-i
         endif !
        enddo !i
!print*,iidt
    do i=1,ipar
         iid(i)=1
         open(1,file=cpathin(1:len_trim(cpathin))//'rdy_clim_71_00.grd',form='binary')
         if(iidt(i).eq.1) then !first day
          do while(.true.)
           read(1)avstid(i,iid(i)),avlat(i,iid(i)),avlon(i,iid(i)),avtime(i,iid(i)),avnlev(i,iid(i)),avnflag(i,iid(i))
           if(avnlev(i,iid(i)).eq.0) goto 333
           read(1)avrf(i,iid(i))  
           iid(i)=iid(i)+1
          enddo !while
          333 continue
         else !not the first day
          iavlev=0
          do while(.true.)
           read(1)avstid(i,iid(i)),avlat(i,iid(i)),avlon(i,iid(i)),avtime(i,iid(i)),avnlev(i,iid(i)),avnflag(i,iid(i))
            ! print*,avstid(i,iid(i)),avlat(i,iid(i)),avlon(i,iid(i)),avtime(i,iid(i)),avnlev(i,iid(i)),avnflag(i,iid(i))
           if(avnlev(i,iid(i)).eq.0) then
            iavlev=iavlev+1
           else
            read(1)avrf(i,iid(i))
           endif !avnlev
            !print*,iavlev
           if(iavlev.eq.iidt(i)-1)  goto 444          
          enddo !while
      444 continue
!          print*,iavlev
      do while(.true.)
           read(1)avstid(i,iid(i)),avlat(i,iid(i)),avlon(i,iid(i)),avtime(i,iid(i)),avnlev(i,iid(i)),avnflag(i,iid(i))
         !  print*,avstid(i,iid(i)),avlat(i,iid(i)),avlon(i,iid(i)),avtime(i,iid(i)),avnlev(i,iid(i)),avnflag(i,iid(i))
           if(avnlev(i,iid(i)).eq.0) then
            goto 555
           else
            read(1)avrf(i,iid(i))
            iid(i)=iid(i)+1
           endif !avnlev
          enddo !while
      555 continue
         endif !iidt(i)
         close(1)
        enddo !i

        irfav=0
        do i=1,ipar
         do j=1,nstation(31)
          do k=1,iid(i)
!           print*,avstid(i,k)             
           write(stid,'(A1,A5,A2)')' ',name(j,31)(1:len_trim(name(j,31))),'  '
          
!          print*,avstidd
       
!          print*,len_trim(avstid(i,k)),len_trim(stid)
!          write(*,'(A)')stid
!          write(*,'(A)')avstid(i,k)
        !   print*,avstid(i,k)
        !   print*,name(j,31)

!           print*,avstid(i,k),stid
!
!          print*,avstid(i,k).eq.stid.and.(avrf(i,k).ne.99999.00)
           if(avstid(i,k).eq.stid.and.(avrf(i,k).ne.99999.00)) then
        !    print*,'eq'
        !        print*,rfrav(j)
            rfrav(j)=rfrav(j)+avrf(i,k)
        !        print*,rfrav(j)

                irfav(j)=irfav(j)+1
           endif !avstid
          enddo !k
         ! print*,name(j,31)
         enddo !j
        enddo !i
!        print*,nstation(31)

        do i=1,nstation(31)
         if(irfav(i).ne.0) then
!          rfrav(i)=rfrav(i)/irfav(i)
!          print*,rfr(i)
      if(rfr(i).ne.-999.0) rfr(i)=rfr(i)-rfrav(i)         
!          print*,rfr(i)
         else
          rfr(i)=-999.0
!          print*,rfr(i)
         endif !irfav
        enddo !

    !=====================================================       
    !-------------------
    !write the data
    open(2,file=cpathout(1:len_trim(cpathout))//ctime(31)//'_'//cpar//'_PT.grd',form='binary') !Form ='unformatted',recordtype='stream')
        do i=1,nstation(31)
         write(stid,'(A5,A3)')name(i,31),'   '
!         print*,stid,name(i,31)
         tim=0.0
         nlev=1
         nflag=1
         write(2)stid,lat(i,31),lon(i,31),tim,nlev,nflag
         write(2)rfr(i)
        enddo !i
        nlev=0
        write(2)stid,lat(nstation(31),31),lon(nstation(31),31),tim,nlev,nflag
    close(2)
        !-------------------
        !write the ctl
        open(3,file=cpathout(1:len_trim(cpathout))//ctime(31)//'_'//cpar//'_PT.ctl')
        write(3,'(A)')'dset '//cpathout(1:len_trim(cpathout))//ctime(31)//'_'//cpar//'_PT.grd'
        write(3,'(A)')'dtype station'
        write(3,'(A)')'stnmap '//cpathout(1:len_trim(cpathout))//ctime(31)//'_'//cpar//'_PTrain.map'
        write(3,'(A)')'title avrage_'//cpar//'_dy rainfall'
        write(3,'(A)')'undef -999.0'
        write(3,'(A)')'tdef 1 linear 18z19Aug2009 1dy'
    write(3,'(A)')'vars 1'
        write(3,'(A)')'rf 0  99 rainfall'
        write(3,'(A)')'endvars'
        close(3)   
!-----------------------------------------------------------------------
!creat gs
open(1,file=cpathriver(1:len_trim(cpathriver))//ctime(31)//'_'//cpar//'_PT.gs')
write(1,'(A)')'"!stnmap -i '//cpathout(1:len_trim(cpathout))//ctime(31)//'_'//cpar//'_PT.ctl"' !lu jing she zhi cheng fan xie gang
write(1,'(A)')'"reinit"'
write(1,'(A)')'"open '//cpathriver(1:len_trim(cpathriver))//'china.ctl"'
write(1,'(A)')'"open '//cpathout(1:len_trim(cpathout))//ctime(31)//'_'//cpar//'_PT.ctl"'
write(1,'(A)')'"set mpdset mres"'
write(1,'(A)')'"set grid off"'
write(1,'(A)')'"set grads off"'
write(1,'(A)')'"set timelab off"'
write(1,'(A)')'"set xlopts 1 4 0.15"'
write(1,'(A)')'"set ylopts 1 4 0.15"'
write(1,'(A)')'"set lat 0 60"'
write(1,'(A)')'"set lon 60 180"'
write(1,'(A)')'"set gxout shaded"'
write(1,'(A)')'"set clevs 0 10 25 50 100"'
write(1,'(A)')'"set ccols 0 10 12 8 7 2"'
write(1,'(A)')'"set csmooth on"'
write(1,'(A)')'"Define A=oacres(g,rf.2,10,5,2,1)"'
write(1,'(A)')'"Define B=smth9(A)"'
!write(1,51)'"Define AA=maskout(A,g.2(t=1)-0.5)"'
!51 format(A35)
write(1,'(A)')'"d B"'
write(1,'(A)')'"cbar"'
write(1,'(A)')'"set map 2 1 6"'
write(1,'(A)')'"draw map"'
write(1,'(A)')'"run '//cpathriver(1:len_trim(cpathriver))//'river_2.gs"'
write(1,'(A)')'"q time"'
write(1,'(A)')'"draw title '//ctime(31)//'_'//cpar//'_PT"'
write(1,'(A)')'"printim '//cpath_pic(1:len_trim(cpath_pic))//ctime(31)//'_'//cpar//'_PT.gif white"'
write(1,'(A)')'"close 2"'
write(1,'(A)')'"close 1"'
close(1)
!-----------------------------------------------------------

   endif !position(31)
  enddo !iday
enddo !imonth
enddo !iyear

end !program





   
   


  

kongfeng0824

自己先做个沙发

mofangbao

介绍一下数据源和大致思路吧 这样别人不仔细看程序就可以大概判断是否需要~

kongfeng0824

mofangbao 发表于 2013-6-7 09:35

                               
登录/注册后可看大图

介绍一下数据源和大致思路吧 这样别人不仔细看程序就可以大概判断是否需要~

数据源就是中国气象数据共享的站点数据。我是算的降水。
然后用grads画图。
其实程序没有什么难的算法，就是简单的批量读取、处理、绘图。一气呵成。
这里的批量读取和清风的批处理的不一样，没有调用cmd。
而且这里面有润平年的区分和大小月的区分。
只要是用过fortran的人稍加看一下就明白。
一段一段的分开来了。

宋秀鹏

了解一下数据源

空谷幽竹

好东西 顶一个

fairyes

请问下，这个批处理文件可以直接在opengrads里运行吗？还是必须要经过fortran compile 软件编译后才能使用？我笔记本上win7系统没法安装编译软件，这时候该怎么办呢？谢谢~

kongfeng0824

fairyes 发表于 2013-6-7 11:18

                               
登录/注册后可看大图

请问下，这个批处理文件可以直接在opengrads里运行吗？还是必须要经过fortran compile 软件编译后才能使用？ ...

fortran的程序

zihan6688

先给个说明，大家好把握需求

liutaoz21

不错，感谢分享。有详细说明的话更好~