短期气候实习报告六

更新时间:2024-03-18 01:29:01 阅读量: 综合文库 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

南京信息工程大学短期气候实验(实习)报告

实验名称 夏季区域降水的定量预测 日期 2017.6.6得分 指导教师 系 大气科学 专业 大气科学 年级 2014级 班次 2姓名 车楚玉 学号 20141301043

一、目的要求:

目的 :掌握短期气候预测中物理统计预测的基本步骤。

要求 :能运用提供的资料和方法子程序,编写或补充完成程序当中的部分片断,了解区域降水的预测方法及其建立过程,输出实验要求的相应结果,并就方法对区域降水的拟合及试验预测效果进行分析 。 二、资料和方法: 资料:

1、前期1月的Nino3.4指数(来自CPC);

2、西太平洋副高脊线、西太平洋副高西伸脊点、亚洲极涡面积、南方涛动指数(来自中国气象局整编的74个环流指数);

3、夏季华北区域10站的降水量,距平百分率。 方法:

回归分析(mregrssion.for)是用来寻找若干变量之间统计关系的一种方法,利用所找到的统计关系对某一变量作出未来时刻的估计,称为回归预报值。 三、实习步骤:

① 说明所用资料方法;

② 计算方法的简单介绍;

③ 输出反映回归效果的参数及回归系数,并就相关参数分析回归效果;

④预测量与回归方程计算的估计值和观测值的历年曲线变化图(1952~2001年),并附简单的说明;

⑤输出独立预测试验的观测与预测值。

四、结果(图表并解释说明):

程序如下:

PROGRAM MAIN

INTEGER,PARAMETER::N=50 INTEGER,PARAMETER::K=5 REAL,DIMENSION(K,N)::X REAL,DIMENSION(N)::Y REAL,DIMENSION(K+1)::A REAL,DIMENSION(K+1,K+1)::B REAL,DIMENSION(K)::V

REAL Q,S,R,U,ind(6,60),year(50),y1(50),pre(7) integer i,j

C x(5,50) y(50) a(6) b(6,6),v(5) C OPEN THE INPUT DATA FILE

OPEN(10,FILE='D:\\duanqi\\shixi6\\shixi.txt')

C READ THE DATA and give data to X and Y

MM=K+1

call DYHG(X,Y,K,MM,N,A,Q,S,R,V,U,B,DYY) write(*,88) A(1) do 89 j=2,MM do i=1,50

read(10,*) year(i),x(1:5,i),y(i) enddo

88 format(/1x,'b 0=',f19.5) 89 write(*,100) j-1,A(j) 100 format(1x,'b',i2,'=',f9.5)

open(11,FILE='D:\\duanqi\\shixi6\\guji.txt')

open(12,FILE='D:\\duanqi\\shixi6\\guji.grd',form='binary') open(13,FILE='D:\\duanqi\\shixi6\\prediction.txt')

do j=1,50

y1(j)=a(1)+a(2)*x(1,j)+a(3)*x(2,j)+a(4)*x(3,j)+a(5)*x(4,j)+a(6) * *x(5,j)

enddo

pre(1)=a(1)+a(2)*26.50+a(3)*13.00+a(4)*100.00+a(5)*196.00+a(6)* * 2.00

pre(2)=a(1)+a(2)*27.76+a(3)*12.00+a(4)*90.0+a(5)*197.00+a(6)*(-1.00) pre(3)=a(1)+a(2)*26.74+a(3)*12.0+a(4)*120.00+a(5)*199.00+a(6)*

*(-11.00)

pre(4)=a(1)+a(2)*27.10+a(3)*13.00+a(4)*110.00+a(5)*220.00+a(6)*3.00 pre(5)=a(1)+a(2)*25.64+a(3)*14.00+a(4)*105.00+a(5)*215.00+a(6)*13.00 pre(6)=a(1)+a(2)*27.26+a(3)*15.00+a(4)*90.00+a(5)*185.00+a(6)* pre(7)=a(1)+a(2)*24.71+a(3)*17.00+a(4)*125.00+a(5)*222.00+a(6)*13.00

print*,pre

*(-8.00)

write(11,*)(y(i),y1(i),i=1,50)

write(12)(y(i),y1(i),i=1,50)

write(13,*)(pre(i),i=1,7)

cccccccccccccccccccccccccccccccccccccccccccccccccccccc

write(*,20)Q,S,R write(*,22)U,DYY

write(*,30)(i,V(i),i=1,K) write(*,40)U

open(6,file='table') ! output data write(6,180) write(6,88) A(1) do 189 j=2,MM write(6,200)

20 format(1x,'Q=',f13.6,3x,'S=',f13.6,3x,'R=',f13.6) 22 format(1x,'U=',f13.6,3x,'DYY=',f13.6) 30 format(1x,'V(',i2,')=',f13.6) 40 format(1x,'U=',f13.6)

180 format(/2x,'regression coefficients:')

189 write(6,100) j-1,A(j)

200 format(/1x,'Generic Analysis of Variance Table for the Multiple * Linear Regression')

write(6,202)

202 format(/1x,'----------------------------------------------------- *---------------')

write(6,204) write(6,202) write(6,206) u2=U/real(K) write(6,208)

K,U,U2 N-1,DYY

204 format(/3x,'Source df SS MS')

206 format(/1x,'Total n-1=',i2,' SST=',f13.4)

208 format(/1x,'Regression K=',i2,' SSR=',f13.4,' MSR=SSR/K=' *,f13.4)

q2=q/real(n-k-1) write(6,209)

n-k-1,q,q2

209 format(/1x,'Residual n-k-1=',i2,' SSE=',f13.4,' MSE=SSE/(n-k-1) *=',f13.4)

f=(U/real(K))/(Q/real(N-K-1)) write(6,220) f write(6,202) close(6) stop end

dimension x(m,n),y(n),a(mm),b(mm,mm),v(m) b(1,1)=n do 20 j=2,mm b(1,j)=0.0 do 10 i=1,n b(j,1)=b(1,j) do 50 i=2,mm do 40 j=i,mm b(i,j)=0.0 do 30 k=1,n b(j,i)=b(i,j)

220 format(/1x,' F=MSR/MSE=',f13.4)

subroutine DYHG(x,y,m,mm,n,a,q,s,r,v,u,b,dyy)

10 b(1,j)=b(1,j)+x(j-1,i) 20 continue

30 b(i,j)=b(i,j)+x(i-1,k)*x(j-1,k) 40 continue 50 continue

a(1)=0.0 do 60 i=1,n do 80 i=2,mm a(i)=0.0 do 70 j=1,n

60 a(1)=a(1)+y(i)

70 a(I)=a(i)+x(i-1,j)*y(j) 80 continue

call CHOLESKY(b,mm,1,a,l) yy=0.0 do 90 i=1,n q=0.0 dyy=0.0 u=0.0 do 110 i=1,n p=a(1)

90 yy=yy+y(i)/n

cccccccccccccccccccccccccccccccccc

do 100 j=1,m

q=q+(y(i)-p)*(y(i)-p) dyy=dyy+(y(i)-yy)*(y(i)-yy) u=u+(yy-p)*(yy-p)

100 p=p+a(j+1)*x(j,i)

110 continue

cccccccccccccccccccccccccccccccccc

subroutine CHOLESKY(a,n,m,d,l) ! Perform the CHOLESKY Decomposition dimension a(n,n),d(n,m) l=1

if(a(1,1)+1.0.eq.1.0)then l=0

write(*,30) return endif

a(1,1)=sqrt(a(1,1)) do 10 j=2,n do 100 i=2,n do 20 j=2,i

if(a(i,i)+1.0.eq.1.0)then l=0

write(*,30) return s=sqrt(q/n) r=sqrt(1.0-q/dyy) do 150 j=1,m p=0.0 do 140 i=1,n pp=a(1) do 130 k=1,m

if(k.ne.j)pp=pp+a(k+1)*x(k,i) p=p+(y(i)-pp)*(y(i)-pp) v(j)=sqrt(1.0-q/p) return end

130 continue 140 continue 150 continue

10 a(1,j)=a(1,j)/a(1,1)

20 a(i,i)=a(i,i)-a(j-1,i)*a(j-1,i)

本文来源:https://www.bwwdw.com/article/eqf8.html

Top