第二章诊断分析1.docx
- 文档编号:26422907
- 上传时间:2023-06-19
- 格式:DOCX
- 页数:5
- 大小:17.46KB
第二章诊断分析1.docx
《第二章诊断分析1.docx》由会员分享,可在线阅读,更多相关《第二章诊断分析1.docx(5页珍藏版)》请在冰豆网上搜索。
第二章诊断分析1
第二章诊断分析1
第二章诊断分析 3.2散度和涡度计算程序 功能 计算散度和涡度 方法说明子程序语句 计算散度的子程序语句 SUBROUTINEDIVER(U,V,P,DL,DM,F0,LL,MM,NN,DIV)SUBROUTINEVORTICITY(U,V,DL,DM,F0,LL,MM,NN,VOR) U——输入参量,实型三维数组,大小为,风速的东西分量,单位:
米/秒。
V——输入参量,实型三维数组,大小为,风速的南北分量,单位:
米/秒。
P——输入参数,实型一维数组,大小为NN,各层的气压值,单位:
百帕。
DL——输入参数,实型常数,X方向的格距。
DM——输入参数,实型常数,Y方向的格距。
F0——输入参数,实型常数,为J=1处的纬度值。
LL——输入参数,整型常数,东西方向的格点数。
MM——输入参数,整型常数,南北方向的格点数。
NN——输入参数,整型常数,垂直方向的层数。
DIV——输出参数,实型三维数组,大小为,散度值,单位:
秒-1。
VOR——输出参数,实型三维数组,大小为,涡度值,单位:
秒-1。
计算散度的子程序语句哑元说明 子程序 计算散度子程序 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,ESPREAL(8),DIMENSION(N):
:
P !
地球半径 REAL(8):
:
DL,DM,F0REAL(8),PARAMETER:
:
AA=!
计算未订正的散度L1=L-1 M1=M-1DOK=1,NDOJ=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 ENDDOENDDO CALLBOUND(DIV,L,M,N)!
进行散度订正 DOJ=1,MDOI=1,LDV(I,J)=DVV(I,J)=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.ENDDOENDDOENDDO 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))ENDDOENDDOENDDORETURNEND SUBROUTINEBOUND(A,L,M,N)REAL(8),DIMENSION(L,M,N):
:
AL1=L-1M1=M-1DOK=1,NDOI=2,L1 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))* A(I,1,K)=2*A(I,2,K)-A(I,3,K) A(L,1,K)=A(L,2,K)+A(L-1,1,K)-(A(L,3,K)+A(L-2,1,K))*A(1,M,K)=A(1,M-1,K)+A(2,M,K)-(A(1,M-2,K)+A(3,M,K))* A(L,M,K)=A(L,M-1,K)+A(L-1,M,K)-(A(L,M-2,K)+A(L-2,M,K))*ENDDOEND SUBROUTINEVORTICITY(U,V,DL,DM,F0,L,M,N,VOR)REAL(8),DIMENSION(L,M,N):
:
U,VREAL(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= REAL(8):
:
DL,DM,F0 L1=L-1M1=M-1DOK=1,NDOJ=2,M1DOI=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)/AAENDDOENDDOENDDO !
地球半径 CALLBOUND(VOR,L,M,N) RETURN END 以1988年5月1日60°E—180°—160°W,0°—70°N范围内的u、v场,计算散度和涡 计算散度和涡度的例子 度,资料共7层,网格点为°×°。
PROGRAMDIVVOR INTEGER,PARAMETER:
:
LL=57INTEGER,PARAMETER:
:
MM=29 INTEGER,PARAMETER:
:
NN=7 REAL(8),DIMENSION(LL,MM,NN):
:
U,VREAL(8),DIMENSION(LL,MM,NN):
:
DIV,VORREAL(8),DIMENSION(NN):
:
P REAL(8),PARAMETER:
:
PI=REAL(8):
:
DL=*PI/180REAL(8):
:
DM=*PI/180 REAL(8):
:
F0=0*PI/180 DATAP/1000,850,700,500,300,200,100/L0=LL OPEN(8,FILE=‘‘) READ(8,’()’)(((U(I,J,K),I=1,LL),J=1,MM),K=1,NN)CLOSE(8) OPEN(9,FILE=‘‘) READ(9,’()’)(((V(I,J,K),I=1,LL),J=1,MM),K=1,NN)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=‘‘)DOK=1,NN WRITE(12,’()’)((DIV(I,J,K),I=1,LL),J=1,MM)ENDDO CLOSE(12) OPEN(13,FILE=‘‘) DOK=1,NN WRITE(13,’()’)((VOR(I,J,K),I=1,LL),J=1,MM)ENDDOCLOSE(13)END 计算结果见图和图:
这里用的都是双精度,即REAL(8),如只用单精度,将其改为REAL(4)即可。
附注 图500hPa散度场 -5-1 图500hPa涡度场 3.3垂直速度ω的计算程序 功能 用运动学方法、准地转ω方程和多层非线性ω平衡方程,求解垂直运动速度。
方法说明 本程序包括两种方法求解垂直运动速度。
运动学方法准地转ω方程子程序语句 运动学法求OMEGA YUNDF(U,V,P,DL,DM,F0,L,M,N,OMEGA)求解准地转OMEGA方程 OME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA) 哑元说明 L——输入参量,整型变量,东西方向格点数。
M——输入参量,整型变量,南北方向格点数。
N——输入参量,整型变量,铅垂方向层数。
U——输入参量,实型三维数组,大小为,风速的东西分量,单位:
米/秒。
V——输入参量,实型三维数组,大小为,风速的南北分量,单位:
米/秒。
T——输入参量,实型三维数组,大小为,温度,单位:
℃。
P——输入参量,实型一维数组,大小为N,铅垂方向N层的气压值,单位:
百帕。
F0——输入参量,实型变量,南边界的纬度。
DL——输入参量,实型变量,东西方向格点间距。
DM——输入参量,实型变量,南北方向格点间距。
OMEGA——输出参量,实型三维数组,大小为,计算出的OMEGA值,单位:
百帕/秒。
子程序 运动学法求OMEGA的子程序 SUBROUTINEYUNDF(U,V,P,DL,DM,F0,L,M,N,OMEGA)REAL(8),DIMENSION(L,M,N):
:
U,V,DIV,OMEGA,OMEGA1REAL(8),DIMENSION(L,M,N):
:
OMTREAL(8),DIMENSION(N):
:
P !
OMT为热力学法算出的顶层OMEGA值 REAL(8):
:
F0,DL,DM OMT=0!
这里取OMT=0,如有资料可用热力学法计算OMEGA1=0 CALLDIVER(U,V,P,DL,DM,F0,L,M,N,DIV)DOI=1,LDOJ=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 ENDDOENDDOEND 求解准地转OMEGA方程的子程序 SUBROUTINEOME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA)!
求解准地转OMEGA方程REAL(8),DIMENSION(L,M,N):
:
U,V,FAI,KAP,SIGMAREAL(8),DIMENSION(L,M,5):
:
T,TP REAL(8),DIMENSION(L,M,N):
:
OMEGA1,OMEGA2,OMEGA,OMEGA0REAL(8),DIMENSION(L,M,N):
:
QA1,QA2,QA3,QB1,QB2,QB3REAL(8),DIMENSION(N):
:
P,SIGMAVREAL(8),DIMENSION(M):
:
FCO REAL(8):
:
EPS,CP,RD,F0,DL,DM,ALFEPS=!
迭代要求的精度ALF=RD=287CP=1005OMEGA1=0 OMEGA2=0 CALLLHS(U,V,P,L,M,N,F0,DL,DM,FAI,KAP)DOJ=1,M!
科氏参数 FCO(J)=2*()*SIN(F0+(J-1)*DM)ENDDO !
中间变量 !
计算静力稳定度 CALLAP1(T,L,M,N,P,TP)DOK=1,NDOJ=1,M DOI=1,L SIGMA(I,J,K)=-(RD/P(K)*(TP(I,J,K)-RD/CP*T(I,J,K)/P(K))) ENDDOENDDOENDDODOK=1,NDOI=1,L DOJ=1,MSIGMAV(K)=SIGMAV(K)+SIGMA(I,J,K) SIGMAV(K)=0 ENDDO ENDDO ENDDODOK=1,N DOJ=1,M DOI=1,L IF(SIGMA(I,J,K) SIGMA(I,J,K)=SIGMAV(K) ENDIF ENDDO ENDDO!
ENDDO 计算右边第一项 CALLLAPLACE(FAI,L,M,N,F0,DL,DM,QA1)DOK=1,NDOJ=1,M DOI=1,L QA1(I,J,K)=QA1(I,J,K)+FCO(J) ENDDOENDDO ENDDO CALLJACOB(FAI,QA1,DL,DM,F0,L,M,N,QA2)CALLAP1(QA2,L,M,N,P,QA3)QA3=QA3*FCO(M/2)N0=0DO 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=‘‘) WRITE(2,’()’)(((OMEGA1(I,J,K),I=1,L),J=1,M),K=1,N)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)*QB3N1=0DO N1=N1+1 OMEGA0=OMEGA2 CALLERROR(OMEGA2,QB3,P,SIGMA,L,M,N,F0,FCO,DL,DM,ALF)IF(MAXVAL(ABS(OMEGA2-OMEGA0)) open(3,file=‘‘) WRITE(3,’()’)(((OMEGA2(I,J,K),I=1,L),J=1,M),K=1,N)close(3) OMEGA=OMEGA1+OMEGA2END SUBROUTINEAP1(A,L,M,N,P,FF)REAL(8),DIMENSION(L,M,N):
:
AREAL(8),DIMENSION(L,M,N):
:
FFREAL(8),DIMENSION(N):
:
PDOJ=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,QAREAL(8):
:
F0,DL,DM,AAAA=L1=L-1M1=M-1DOK=1,NDOI=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 ENDDOENDDO 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,JCREAL(8):
:
F0,DL,DM,J1,J2 REAL(8):
:
AA=L1=L-1M1=M-1DOK=1,NDOJ=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)=*J1+*J2ENDDO ENDDOENDDO 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,SIREAL(8),DIMENSION(N):
:
PREAL(8),DIMENSION(M):
:
FCOREAL(8):
:
F0,DL,DM,AA,EPS,ALF,REPS=AA=L1=L-1M1=M-1 N1=N-1 DOK=2,N1DOI=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)ENDDOENDDOENDDOEND
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第二 诊断 分析