请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 6083|回复: 10

[脚本编辑] 紧急求助!画160站降水距平图

[复制链接]

新浪微博达人勋

发表于 2018-6-1 21:05:47 | 显示全部楼层 |阅读模式

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

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

x
论文需要画这个图,我是做的平流层QBO对我国夏季降水的影响,找到了洪hong的一个帖子:中国160站降水的画图
http://bbs.06climate.com/forum.php?mod=viewthread&tid=44713&fromuid=97660
(出处: 气象家园)
楼主给的程序包很有用写的很好,图做出来真的很漂亮,但是楼主的程Fortran序是做的指数与降水距平的相关系数,我只需要做降水距平图,不知道该怎么改,这个程序算出的降水距平是二维数组,楼主给的绘图程序是三个步骤,第一个计算,第二步会生成一个grb文件,第三部gs直接绘图,但是想下面是我改了一半的第一步程序,想算1979-2017年的西风位相年份的降水距平,请各位大佬救救命

program rain_160_relation
integer,parameter :: m=160,nt=67,n=39,wt=21   !1951年到2017年=67年 我用的是1979年到2017年=39年 1979年=第29格点开始
integer it
real a6(m,nt),a7(m,nt),a8(m,nt),w6(m,wt),w7(m,wt),w8(m,wt),ra(m,wt),rsum(m),rain(m,wt)
open(1,file='C:\data\zother\fortran_160_1\r1606.txt')
do it=1,nt
read(1,*)(a6(i,it),i=1,m)
enddo

open(2,file='C:\data\zother\fortran_160_1\r1607.txt')
do it=1,nt
read(2,*)(a7(i,it),i=1,m)
enddo

open(3,file='C:\data\zother\fortran_160_1\r1608.txt')
do it=1,nt
read(3,*)(a8(i,it),i=1,m)
enddo

open(4,file='C:\data\zother\fortran_160_1\ww6.txt')
do it=1,wt
read(4,*)(w6(i,it),i=1,m)
enddo

open(5,file='C:\data\zother\fortran_160_1\ww7.txt')
do it=1,wt
read(5,*)(w7(i,it),i=1,m)
enddo

open(6,file='C:\data\zother\fortran_160_1\ww8.txt')
do it=1,wt
read(6,*)(w8(i,it),i=1,m)
enddo






!------------------------------------求降水距平百分率-ra(m,n)----------------------------------------------
do i=1,160
      rsum(i)=0.0
        do it=29,67
        rsum(i)=rsum(i)+a6(i,it)+a7(i,it)+a8(i,it)
        enddo

        rsum(i)=rsum(i)/39.0
        do it=1,21
        ra(i)=w6(i,it)+w7(i,it)+w8(i,it)-rsum(i))*100/rsum(i)
        enddo
        enddo
!------------------------------------求降水距平场-rain(m,n)----------------------------------------------
do i=1,160
      rsum(i)=0.0
        do it=29,67
        rsum(i)=rsum(i)+a6(i,it)+a7(i,it)+a8(i,it)
        enddo

        rsum(i)=rsum(i)/39.0
        do it=1,21
        rain(i,it)=w6(i,it)+w7(i,it)+w8(i,it)-rsum(i)
        enddo
        enddo








!call relation(ra,SMI,rrSMI,m,n)


open(8,file='rrSMI.txt')
do it=1,m
write(8,*)()
enddo



end
!--------------------------------------求相关系数子程序-----------------------------------------------
subroutine relation(rain,heat,r,m0,n0)
integer m0,n0
real rain(m0,n0),heat(n0),r(m0)
real avr(160),avh,sumh,sumr(160),sumk(160),suml(160),sump(160)
avr=0.0
avh=0.0
sumh=0.0
sumr=0.0
sumk=0.0
suml=0.0
sump=0.0
        do i=1,n0
                sumh=sumh+heat(i)
                enddo
                avh=sumh/n0
           do i=1,m0
           do j=1,n0
           sumr(i)= sumr(i)+rain(i,j)
           enddo
           avr(i)=sumr(i)/n0
           enddo
           
           do i=1,m0
           do j=1,n0
          sumk(i)=sumk(i)+(rain(i,j)-avr(i))*(heat(j)-avh)
          suml(i)=suml(i)+(rain(i,j)-avr(i))**2
          sump(i)=sump(i)+(heat(j)-avh)**2
          enddo
          r(i)=(sumk(i)/n0)/((sqrt(suml(i))*sqrt(sump(i)))/n0)
          enddo
end subroutine

这是第二步程序
program ueddcu
real lat(160),lon(160),rrSIM(160)
character*8 id(160)        


open(1,file='C:\data\zother\fortran_zhandian\lat_lon.txt')
        do i=1,160
        read(1,*)lat(i),lon(i)
        enddo

        open(2,file='C:\data\zother\fortran_zhandian\rrSMI.txt')        !只修改文件名
        do i=1,160
        read(2,*)rrsim(i)
        enddo

open(8,file='C:\data\zother\fortran_zhandian\rrSIM.grb',form='binary')

!ccccccccccccccccccc写站点数据

        do j=1,160
        id(j)=char(j)
        tim=0.0
        nlev=1
        nflag=1
        write(8)id(j),lat(j),lon(j),tim,nlev,nflag,rrsim(j)
        enddo
        tim=0.0
        nlev=0
        nflag=1
        write(8)id(j-1),lat(j-1),lon(j-1),tim,nlev,nflag


        end
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-6-1 21:12:25 | 显示全部楼层
是某一年的降水百分率?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-6-1 21:14:21 | 显示全部楼层
本帖最后由 高小一 于 2018-6-1 21:17 编辑
qxsw2016 发表于 2018-6-1 21:12
是某一年的降水百分率?

是1979-2017年的,,那个r1606,r1607,r1608是1951到2017的,ww6,ww7,ww8是79年至17年间西风位相年份的降水数据,想做39年一起的降水距平,还有QBO东风位相和西风位相年份的降水距平,这个贴出来的打算算西风位相年份的降水距平的
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-6-1 21:21:26 | 显示全部楼层
ww6、7、8存的是什么
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-6-1 21:24:19 | 显示全部楼层
是79年至17年间西风位相年份的降水数据,6,7,8月的月平均
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-6-1 21:36:15 | 显示全部楼层
数据都在这39年里吧?第一步试试下面代码。westphase(21) 把西位相年赋初值
program rain_160_relation
integer,parameter :: m=160,nt=67,n=39,wt=21   !1951年到2017年=67年 我用的是1979年到2017年=39年 1979年=第29格点开始
integer it,i,j
integer westphase(21) !比如东位像年数组
real a6(m,nt),a7(m,nt),a8(m,nt),w6(m,wt),w7(m,wt),w8(m,wt),ra(m,wt),rsum(m),rain(m,wt)
!应该是1951-2017年6、7、8月降水
open(1,file='C:\data\zother\fortran_160_1\r1606.txt')
do it=1,nt
read(1,*)(a6(i,it),i=1,m)
enddo

open(2,file='C:\data\zother\fortran_160_1\r1607.txt')
do it=1,nt
read(2,*)(a7(i,it),i=1,m)
enddo

open(3,file='C:\data\zother\fortran_160_1\r1608.txt')
do it=1,nt
read(3,*)(a8(i,it),i=1,m)
enddo


do i=1,160
      rsum(i)=0.0
        do it=29,67
        rsum(i)=rsum(i)+a6(i,it)+a7(i,it)+a8(i,it)
        enddo

        rsum(i)=rsum(i)/39.0  !39年平均值
enddo


do i= 1,m
   do j = 1,wt
     it = westphase(j)-1950
     ra(i) = ra(i)+a6(i,it)+a7(i,it)+a8(i,it)-rsum(i)
  enddo
    ra(i) = ra(i)/wt
enddo
end

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-6-1 22:59:53 | 显示全部楼层
本帖最后由 高小一 于 2018-6-1 23:05 编辑

我把计算的部分改成这样了,然后write ra(m),
do i=1,160
      rsum(i)=0.0
        do it=29,67
        rsum(i)=rsum(i)+a6(i,it)+a7(i,it)+a8(i,it)
        enddo

        rsum(i)=rsum(i)/39.0
        enddo
        ra(i)=0
        do i=1,m
          do j=1,wt

        ra(i)=ra(i)+w6(i,j)+w7(i,j)+w8(i,j)-rsum(i)
        enddo
        ra(i)=ra(i)/21
        
        enddo
open(8,file='xfjp.txt')
do it=1,m
write(8,*)ra(it)
enddo

      

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-6-1 23:11:43 | 显示全部楼层
如果要算降水距平百分率要怎么改啊
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-6-1 23:33:59 | 显示全部楼层
想算距平百分率,所以又改了改,现在是这样
do i=1,160
      rsum(i)=0.0
        do it=29,67
        rsum(i)=rsum(i)+a6(i,it)+a7(i,it)+a8(i,it)
        enddo

        rsum(i)=rsum(i)/39.0
        enddo
       
        do i=1,m
                ra(i)=0
          do j=1,wt

        ra(i)=ra(i)+(w6(i,j)+w7(i,j)+w8(i,j)-rsum(i))*100/rsum(i)
        enddo

       
        enddo
结果感觉不太对劲,大概是这样
172.726400
       99.611600
       16.171700
       69.484150
      144.060200
       25.788830
       63.967090
      -52.500080
        5.710473
       29.359540
       17.094020
      -34.768710
     -113.244600
     -107.693900
        8.536527
      -20.607120
        7.229273
      -93.274000
     -136.549600
       12.902360
      -38.342630
      -40.204960
      -15.498690
     -132.723900
      -28.912740
       21.477280
     -146.916800
      -77.397480
      -10.772410
       -4.030693
       39.994440
     -225.563900
     -117.609700
      -34.924090
      -47.177360
     -201.077700
     -125.474700
     -215.416800
       17.917530
      -21.492380
     -104.473600
     -213.299100
      -91.894630
      -83.078570
     -119.139400
      -77.495160
     -118.641600
      -88.202130
     -264.271100
     -171.843000
     -213.295500
     -235.622600
      -79.736980
     -188.821500
     -241.488900
       56.026760
       32.184810
     -108.243200
      112.838900
       70.020060
     -139.237900
     -136.642300
       60.934850
        6.321421
       30.765770
        6.430769
       15.457500
       22.019260
       59.014640
      224.060200
      135.881800
       67.185730
      149.038000
      -75.006820
       33.340150
      -31.777480
      105.145100
       39.900060
       66.497090
       38.557220
      -52.947740
      184.382400
      148.812400
      119.391000
       24.117600
      -12.622850
       19.906270
       26.491140
      -30.609750
      -52.361650
       66.260380
      -19.523440
      -36.677890
      -85.490170
       38.482250
     -176.606100
     -120.072200
       39.440190
    2.337913E-01
      -11.226470
       75.857160
      -13.891620
      -67.532550
      -42.891020
      -74.946710
     -171.592800
      -98.786680
       59.705070
      109.490200
     -190.223900
      -48.296760
      -13.481920
       27.763130
      110.776400
       22.376690
       70.226910
      -56.103510
       88.716370
       40.658230
       24.476730
       71.278340
       16.032890
       58.648270
      -92.849630
       52.424240
     -228.624200
       62.074450
       18.450080
       -4.983445
     -208.677400
     -127.510600
      -52.507870
        5.754626
      -11.946910
       10.377240
      177.838300
       78.605670
      -19.801350
       84.175140
     -125.660100
      -90.040690
      -44.690540
      -88.480170
       24.523120
      274.959800
       53.257320
     -172.496500
     -130.616300
     -351.245500
      -32.289370
     -427.123000
     -307.368400
     -531.147600
     -155.718600
     -292.455200
     -207.966700
       15.276780
       90.121820
       84.907520
     -282.822300
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-6-1 23:56:30 | 显示全部楼层
高小一 发表于 2018-6-1 23:33
想算距平百分率,所以又改了改,现在是这样
do i=1,160
      rsum(i)=0.0

ra(i)=ra(i)/21
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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