第二章诊断分析1.docx
- 文档编号:8872691
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:16
- 大小:18.67KB
第二章诊断分析1.docx
《第二章诊断分析1.docx》由会员分享,可在线阅读,更多相关《第二章诊断分析1.docx(16页珍藏版)》请在冰豆网上搜索。
第二章诊断分析1
第二章诊断分析
3.2散度和涡度计算程序
2.2.1功能
计算散度和涡度
2.2.2方法说明
2.2.3子程序语句
2.2.3.1计算散度的子程序语句
SUBROUTINEDIVER(U,V,P,DL,DM,F0,LL,MM,NN,DIV)
2.2.3.2.计算散度的子程序语句
SUBROUTINEVORTICITY(U,V,DL,DM,F0,LL,MM,NN,VOR)
2.2.4哑元说明
U——输入参量,实型三维数组,大小为(LL*MM*NN),风速的东西分量,单位:
米/秒。
V——输入参量,实型三维数组,大小为(LL*MM*NN),风速的南北分量,单位:
米/秒。
P——输入参数,实型一维数组,大小为NN,各层的气压值,单位:
百帕(hPa)。
DL——输入参数,实型常数,X方向的格距(单位为弧度)。
DM——输入参数,实型常数,Y方向的格距(单位为弧度)。
F0——输入参数,实型常数,为J=1处的纬度值(单位为弧度)。
LL——输入参数,整型常数,东西方向的格点数。
MM——输入参数,整型常数,南北方向的格点数。
NN——输入参数,整型常数,垂直方向的层数。
DIV——输出参数,实型三维数组,大小为(LL*MM*NN),散度值,单位:
秒-1。
VOR——输出参数,实型三维数组,大小为(LL*MM*NN),涡度值,单位:
秒-1。
2.2.5子程序
2.2.5.1计算散度子程序(子程序名为:
DIVER)
SUBROUTINEDIVER(U,V,P,DL,DM,F0,L,M,N,DIV)
!
U、V分别为风速的东西和南北分量;P为气压值;DX、DY分别为X、Y方向的格距(单位为弧度);
!
F0为J=1处的纬度值(弧度);L、M分别为X、Y方向的格点数;N为垂直方向的层数;DIV为散度。
REAL(8),DIMENSION(L,M,N):
:
U,V
REAL(8),DIMENSION(L,M,N):
:
DIV
REAL(8),DIMENSION(L,M):
:
DV,DVV,ESP
REAL(8),DIMENSION(N):
:
P
REAL(8):
:
DL,DM,F0
REAL(8),PARAMETER:
:
AA=6.371E6!
地球半径
!
计算未订正的散度
L1=L-1
M1=M-1
DOK=1,N
DOJ=2,M1
DOI=2,L1
DIV(I,J,K)=(U(I+1,J,K)-U(I-1,J,K))/(AA*COS(DM*(J-1)+F0)*2*DL)&
+(V(I,J+1,K)-V(I,J-1,K))/(AA*2*DM)&
-V(I,J,K)*TAN(DM*(J-1)+F0)/AA
ENDDO
ENDDO
ENDDO
CALLBOUND(DIV,L,M,N)
!
进行散度订正
DOJ=1,M
DOI=1,L
DV(I,J)=0.0
DVV(I,J)=0.0
DOK=1,N-1
DV(I,J)=DV(I,J)+(DIV(I,J,K)+DIV(I,J,K+1))*(P(K+1)-P(K))/2.
DVV(I,J)=DVV(I,J)+(ABS(DIV(I,J,K))+ABS(DIV(I,J,K+1)))*(P(K+1)-P(K))/2.
ENDDO
ENDDO
ENDDO
DOJ=1,M
DOI=1,L
ESP(I,J)=-DV(I,J)/DVV(I,J)
DOK=1,N
DIV(I,J,K)=DIV(I,J,K)+ESP(I,J)*ABS(DIV(I,J,K))
ENDDO
ENDDO
ENDDO
RETURN
END
SUBROUTINEBOUND(A,L,M,N)
REAL(8),DIMENSION(L,M,N):
:
A
L1=L-1
M1=M-1
DOK=1,N
DOI=2,L1
A(I,1,K)=2*A(I,2,K)-A(I,3,K)
A(I,M,K)=2*A(I,M-1,K)-A(I,M-2,K)
ENDDO
DOJ=2,M1
A(1,J,K)=2*A(2,J,K)-A(3,J,K)
A(L,J,K)=2*A(L-1,J,K)-A(L-2,J,K)
ENDDO
A(1,1,K)=A(1,2,K)+A(2,1,K)-(A(1,3,K)+A(3,1,K))*0.5
A(L,1,K)=A(L,2,K)+A(L-1,1,K)-(A(L,3,K)+A(L-2,1,K))*0.5
A(1,M,K)=A(1,M-1,K)+A(2,M,K)-(A(1,M-2,K)+A(3,M,K))*0.5
A(L,M,K)=A(L,M-1,K)+A(L-1,M,K)-(A(L,M-2,K)+A(L-2,M,K))*0.5
ENDDO
END
2.2.5.2.计算涡度子程序(子程序名为:
VORTICITY)
SUBROUTINEVORTICITY(U,V,DL,DM,F0,L,M,N,VOR)
REAL(8),DIMENSION(L,M,N):
:
U,V
REAL(8),DIMENSION(L,M,N):
:
VOR
!
U、V分别为风速的东西和南北分量;DX、DY分别为X、Y方向的格距(单位为弧度);
!
F0为J=1处的纬度值(弧度);L、M分别为X、Y方向的格点数;N为垂直方向的层数;
!
VOR为涡度。
REAL(8),PARAMETER:
:
AA=6.371E6!
地球半径
REAL(8):
:
DL,DM,F0
L1=L-1
M1=M-1
DOK=1,N
DOJ=2,M1
DOI=2,L1
VOR(I,J,K)=(V(I+1,J,K)-V(I-1,J,K))/(AA*COS(DM*(J-1)+F0)*2*DL)&
-(U(I,J+1,K)-U(I,J-1,K))/(AA*2*DM)&
+U(I,J,K)*TAN(DM*(J-1)+F0)/AA
ENDDO
ENDDO
ENDDO
CALLBOUND(VOR,L,M,N)
RETURN
END
2.2.6计算散度和涡度的例子
以1988年5月1日60°E—180°—160°W,0°—70°N范围内的u、v场,计算散度和涡度,资料共7层,网格点为2.5°×2.5°。
PROGRAMDIVVOR
INTEGER,PARAMETER:
:
LL=57
INTEGER,PARAMETER:
:
MM=29
INTEGER,PARAMETER:
:
NN=7
REAL(8),DIMENSION(LL,MM,NN):
:
U,V
REAL(8),DIMENSION(LL,MM,NN):
:
DIV,VOR
REAL(8),DIMENSION(NN):
:
P
REAL(8),PARAMETER:
:
PI=3.1415926
REAL(8):
:
DL=2.5*PI/180
REAL(8):
:
DM=2.5*PI/180
REAL(8):
:
F0=0*PI/180
DATAP/1000,850,700,500,300,200,100/
L0=LL
OPEN(8,FILE='U880501.DAT')
READ(8,'(
CLOSE(8)
OPEN(9,FILE='V880501.DAT')
READ(9,'(
CLOSE(9)
CALLDIVER(U,V,P,DL,DM,F0,LL,MM,NN,DIV)
CALLVORTICITY(U,V,DL,DM,F0,LL,MM,NN,VOR)
OPEN(12,FILE='DIV.DAT')
DOK=1,NN
WRITE(12,'(
ENDDO
CLOSE(12)
OPEN(13,FILE='VOR.DAT')
DOK=1,NN
WRITE(13,'(
ENDDO
CLOSE(13)
END
计算结果见图2.2.3和图2.2.3:
(这里给出的图的范围为20°N-50°N,100°E-130°E)
2.2.7附注
这里用的都是双精度,即REAL(8),如只用单精度,将其改为REAL(4)即可。
图2.2.3500hPa散度场(单位:
10-5秒-1)
图2.2.4500hPa涡度场(单位:
10-5秒-1)
3.3垂直速度ω的计算程序
2.3.1功能
用运动学方法、准地转ω方程和多层非线性ω平衡方程,求解垂直运动速度。
2.3.2方法说明
本程序包括两种方法求解垂直运动速度。
2.3.2.1运动学方法
2.3.2.2准地转ω方程
2.3.3子程序语句
2.3.3.1运动学法求OMEGA
YUNDF(U,V,P,DL,DM,F0,L,M,N,OMEGA)
2.3.3.2求解准地转OMEGA方程
OME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA)
2.3.4哑元说明
L——输入参量,整型变量,东西方向格点数。
M——输入参量,整型变量,南北方向格点数。
N——输入参量,整型变量,铅垂方向层数。
U——输入参量,实型三维数组,大小为(L*M*N),风速的东西分量,单位:
米/秒。
V——输入参量,实型三维数组,大小为(L*M*N),风速的南北分量,单位:
米/秒。
T——输入参量,实型三维数组,大小为(L*M*N),温度,单位:
℃。
P——输入参量,实型一维数组,大小为N,铅垂方向N层的气压值,单位:
百帕。
F0——输入参量,实型变量,南边界的纬度。
DL——输入参量,实型变量,东西方向格点间距。
DM——输入参量,实型变量,南北方向格点间距。
OMEGA——输出参量,实型三维数组,大小为(L*M*N),计算出的OMEGA值,单位:
百帕/秒。
2.3.5子程序
2.3.5.1运动学法求OMEGA的子程序
SUBROUTINEYUNDF(U,V,P,DL,DM,F0,L,M,N,OMEGA)
REAL(8),DIMENSION(L,M,N):
:
U,V,DIV,OMEGA,OMEGA1
REAL(8),DIMENSION(L,M,N):
:
OMT!
OMT为由热力学法算出的顶层OMEGA值
REAL(8),DIMENSION(N):
:
P
REAL(8):
:
F0,DL,DM
OMT=0!
这里取OMT=0,如有资料可用热力学法计算
OMEGA1=0
CALLDIVER(U,V,P,DL,DM,F0,L,M,N,DIV)
DOI=1,L
DOJ=1,M
DOK=2,N
OMEGA1(I,J,K)=OMEGA1(I,J,K-1)-(DIV(I,J,K-1)+DIV(I,J,K))/2*(P(K)-P(K-1))
ENDDO
DOK=1,N
OMEGA(I,J,K)=OMEGA1(I,J,K)-K*(K-1)/(N*(N-1))*(OMEGA1(I,J,N)-OMT(I,J,N))
ENDDO
ENDDO
ENDDO
END
2.3.5.2求解准地转OMEGA方程的子程序
SUBROUTINEOME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA)!
求解准地转OMEGA方程
REAL(8),DIMENSION(L,M,N):
:
U,V,FAI,KAP,SIGMA
REAL(8),DIMENSION(L,M,5):
:
T,TP
REAL(8),DIMENSION(L,M,N):
:
OMEGA1,OMEGA2,OMEGA,OMEGA0
REAL(8),DIMENSION(L,M,N):
:
QA1,QA2,QA3,QB1,QB2,QB3!
中间变量
REAL(8),DIMENSION(N):
:
P,SIGMAV
REAL(8),DIMENSION(M):
:
FCO
REAL(8):
:
EPS,CP,RD,F0,DL,DM,ALF
EPS=1.E-5!
迭代要求的精度
ALF=1.44
RD=287
CP=1005
OMEGA1=0
OMEGA2=0
CALLLHS(U,V,P,L,M,N,F0,DL,DM,FAI,KAP)
DOJ=1,M!
科氏参数
FCO(J)=2*(7.29E-5)*SIN(F0+(J-1)*DM)
ENDDO
!
计算静力稳定度
CALLAP1(T,L,M,N,P,TP)
DOK=1,N
DOJ=1,M
DOI=1,L
SIGMA(I,J,K)=-(RD/P(K)*(TP(I,J,K)-RD/CP*T(I,J,K)/P(K)))
ENDDO
ENDDO
ENDDO
DOK=1,N
SIGMAV(K)=0
DOI=1,L
DOJ=1,M
SIGMAV(K)=SIGMAV(K)+SIGMA(I,J,K)
ENDDO
ENDDO
ENDDO
DOK=1,N
DOJ=1,M
DOI=1,L
IF(SIGMA(I,J,K)<=0)THEN
SIGMA(I,J,K)=SIGMAV(K)
ENDIF
ENDDO
ENDDO
ENDDO
!
计算右边第一项
CALLLAPLACE(FAI,L,M,N,F0,DL,DM,QA1)
DOK=1,N
DOJ=1,M
DOI=1,L
QA1(I,J,K)=QA1(I,J,K)+FCO(J)
ENDDO
ENDDO
ENDDO
CALLJACOB(FAI,QA1,DL,DM,F0,L,M,N,QA2)
CALLAP1(QA2,L,M,N,P,QA3)
QA3=QA3*FCO(M/2)
N0=0
DO
N0=N0+1
OMEGA0=OMEGA1
CALLERROR(OMEGA1,QA3,P,SIGMA,L,M,N,F0,FCO,DL,DM,ALF)
IF(MAXVAL(ABS(OMEGA1-OMEGA0)) ENDDO open(2,file='omega1.dat') WRITE(2,'( close (2) ! 计算右边第二项 CALLAP1(FAI,L,M,N,P,QB1) CALLJACOB(FAI,QB1,DL,DM,F0,L,M,N,QB2) CALLLAPLACE(QB2,L,M,N,F0,DL,DM,QB3) QB3=-FCO(M/2)*QB3 N1=0 DO N1=N1+1 OMEGA0=OMEGA2 CALLERROR(OMEGA2,QB3,P,SIGMA,L,M,N,F0,FCO,DL,DM,ALF) IF(MAXVAL(ABS(OMEGA2-OMEGA0)) ENDDO open(3,file='omega2.dat') WRITE(3,'( close(3) OMEGA=OMEGA1+OMEGA2 END SUBROUTINEAP1(A,L,M,N,P,FF)! 计算垂直方向一次不等距差分 REAL(8),DIMENSION(L,M,N): : A REAL(8),DIMENSION(L,M,N): : FF REAL(8),DIMENSION(N): : P DOJ=1,M DOI=1,L DOK=2,N-1 FF(I,J,K)=((P(K)-P(K-1))*A(I,J,K+1)+(P(K+1)+P(K-1)-2*P(K))*A(I,J,K)& -(P(K+1)-P(K))*A(I,J,K-1))/(2*(P(K)-P(K-1))*(P(K+1)-P(K))) ENDDO ENDDO ENDDO END SUBROUTINELAPLACE(A,L,M,N,F0,DL,DM,QA) REAL(8),DIMENSION(L,M,N): : A,QA REAL(8): : F0,DL,DM,AA AA=6.371E6 L1=L-1 M1=M-1 DOK=1,N DOI=2,L1 DOJ=2,M1 QA(I,J,K)=(A(I+1,J,K)+A(I-1,J,K)-2*A(I,J,K))/(DL*COS(F0+(J-1)& *DM))**2+(A(I,J+1,K)+A(I,J-1,K)-2*A(I,J,K))/(DM*DM)& -TAN(F0+(J-1)*DM)*(A(I,J+1,K)-A(I,J-1,K))/(2*DM) ENDDO ENDDO ENDDO QA=QA/(AA*AA) CALLBOUND(QA,L,M,N) END SUBROUTINEJACOB(A,B,DL,DM,F0,L,M,N,JC)! 采用两种JACOB差分方法的平均 REAL(8),DIMENSION(L,M,N): : A,B,JC REAL(8): : F0,DL,DM,J1,J2 REAL(8): : AA=6.371E6 L1=L-1 M1=M-1 DOK=1,N DOJ=2,M1 DOI=2,L1 J1=-((B(I+1,J,K)-B(I-1,J,K))*(A(I,J+1,K)-A(I,J-1,K))& -(B(I,J+1,K)-B(I,J-1,K))*(A(I+1,J,K)-A(I-1,J,K)))& /(AA*AA*COS(F0+(J-1)*DM)*4*DL*DM) J2=-(B(I+1,J+1,K)*(A(I,J+1,K)-A(I+1,J,K))& -B(I-1,J-1,K)*(A(I-1,J,K)-A(I,J-1,K))& -B(I-1,J,K)*(A(I,J+1,K)-A(I-1,J,K))& +B(I+1,J-1,K)*(A(I+1,J,K)-A(I,J-1,K)))& /(AA*AA*COS(F0+(J-1)*DM)*4*DL*DM) JC(I,J,K)=0.6*J1+0.4*J2 ENDDO ENDDO ENDDO CALLBOUND(JC,L,M,N) END SUBROUTINEERROR(A,Q,P,SI,L,M,N,F0,FCO,DL,DM,ALF) REAL(8),DIMENSION(L,M,N): : A,Q,SI REAL(8),DIMENSION(N): : P REAL(8),DIMENSION(M): : FCO REAL(8): : F0,DL,DM,AA,EPS,ALF,R EPS=1.E3 AA=6.371E6 L1=L-1 M1=M-1 N1=N-1 DOK=2,N1 DOI=2,L1 DOJ=2,M1 R=((A(I+1,J,K)+A(I-1,J,K))/(DL*COS(F0+(J-1)*DM))**2& +(A(I,J+1,K)+A(I,J-1,K))/DM**2-TAN(F0+(J-1)*DM)*& (A(I,J+1,K)-A(I,J-1,K))/(2*DM)+(AA*FCO(M/2))**2*& ((P(K)-P(K-1))*A(I,J,K+1)+(P(K+1)-P(K))*A(I,J,K-1))& /(SI(I,J,K)*(P(K)-P(K-1))**2*(P(K+1)-P(K)))& -AA*AA*Q(I,J,K)/SI(I,J,K))/(2/(DL*COS(F0+(J-1)*DM))**2& +2/DM**2+(AA*FCO(M/2))**2*(P(K+1)-P(K-1))& /(SI(I,J,K)*(P(K)-P(K-1))**2*(P(K+1)-P(K)))) A(I,J,K)=ALF*R+(1-ALF)*A(I,J,K) ENDDO ENDDO ENDDO END
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第二 诊断 分析
![提示](https://static.bdocx.com/images/bang_tan.gif)