第二章诊断分析1.docx

资源描述

第二章诊断分析1.docx

《第二章诊断分析1.docx》由会员分享，可在线阅读，更多相关《第二章诊断分析1.docx（16页珍藏版）》请在冰豆网上搜索。

第二章诊断分析1.docx

第二章诊断分析1

第二章诊断分析

３．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）:

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

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

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

RETURN

END

SUBROUTINEBOUND（A,L,M,N）

REAL（8）,DIMENSION（L,M,N）:

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

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）:

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,'（F7.1）'）（（（U（I,J,K）,I=1,LL）,J=1,MM）,K=1,NN）

CLOSE（8）

OPEN（9,FILE='V880501.DAT'）

READ（9,'（F7.1）'）（（（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='DIV.DAT'）

DOK=1,NN

WRITE（12,'（E11.4）'）（（DIV（I,J,K）,I=1,LL）,J=1,MM）

ENDDO

CLOSE（12）

OPEN（13,FILE='VOR.DAT'）

DOK=1,NN

WRITE（13,'（E11.4）'）（（VOR（I,J,K）,I=1,LL）,J=1,MM）

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垂直速度ω的计算程序

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）:

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

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

DOK=1,N

SIGMAV（K）=0

DOI=1,L

DOJ=1,M

SIGMAV（K）=SIGMAV（K）+SIGMA（I,J,K）

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

计算右边第一项

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

CALLJACOB（FAI,QA1,DL,DM,F0,L,M,N,QA2）

CALLAP1（QA2,L,M,N,P,QA3）

QA3=QA3*FCO（M/2）

N0=0

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,'（E11.4）'）（（（OMEGA1（I,J,K）,I=1,L）,J=1,M）,K=1,N）

（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

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,'（E11.4）'）（（（OMEGA2（I,J,K）,I=1,L）,J=1,M）,K=1,N）

close（3）

OMEGA=OMEGA1+OMEGA2

END

SUBROUTINEAP1（A,L,M,N,P,FF）!

计算垂直方向一次不等距差分

REAL（8）,DIMENSION（L,M,N）:

REAL（8）,DIMENSION（N）:

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

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

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

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）:

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

END

展开阅读全文