第二章诊断分析
３．2散度和涡度计算程序
2.2.1功能
计算散度和涡度
2.2.2方法说明 2.2.3子程序语句
2.2.3.1计算散度的子程序语句
SUBROUTINE DIVER(U,V,P,DL,DM,F0,LL,MM,NN,DIV) SUBROUTINE VORTICITY(U,V,DL,DM,F0,LL,MM,NN,VOR)
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.3.2.计算散度的子程序语句 2.2.4哑元说明
2.2.5子程序
2.2.5.1计算散度子程序（子程序名为：DIVER）
SUBROUTINE DIVER(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 DO K=1,N DO J=2,M1
DO I=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 END DO
END DO END DO
CALL BOUND(DIV,L,M,N) ! 进行散度订正
DO J=1,M DO I=1,L DV(I,J)=0.0 DVV(I,J)=0.0 DO K=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. END DO END DO END DO
DO J=1,M DO I=1,L
ESP(I,J)=-DV(I,J)/DVV(I,J) DO K=1,N
DIV(I,J,K)=DIV(I,J,K)+ESP(I,J)*ABS(DIV(I,J,K)) END DO END DO END DO RETURN END
SUBROUTINE BOUND(A,L,M,N) REAL(8),DIMENSION(L,M,N)::A L1=L-1 M1=M-1 DO K=1,N DO I=2,L1
A(I,M,K)=2*A(I,M-1,K)-A(I,M-2,K) END DO
DO J=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) END DO
A(1,1,K)=A(1,2,K)+A(2,1,K)-(A(1,3,K)+A(3,1,K))*0.5
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))*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 END DO END
SUBROUTINE VORTICITY(U,V,DL,DM,F0,L,M,N,VOR) REAL(8),DIMENSION(L,M,N)::U,V REAL(8),DIMENSION(L,M,N)::VOR
2.2.5.2.计算涡度子程序（子程序名为：VORTICITY）
! 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 DO K=1,N DO J=2,M1 DO I=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 END DO END DO END DO
!地球半径
CALL BOUND(VOR,L,M,N) RETURN
END
以1988年5月1日60°E—180°—160°W，0°—70°N范围内的u、v场，计算散度和涡
2.2.6计算散度和涡度的例子
度，资料共7层，网格点为2.5°×2.5°。 PROGRAM DIVVOR
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
DATA P/1000,850,700,500,300,200,100/ L0=LL
OPEN(8,FILE='U880501.DAT')
READ(8,'(
OPEN(9,FILE='V880501.DAT')
READ(9,'(
CALL DIVER(U,V,P,DL,DM,F0,LL,MM,NN,DIV) CALL VORTICITY(U,V,DL,DM,F0,LL,MM,NN,VOR) OPEN(12,FILE='DIV.DAT') DO K=1,NN
WRITE(12,'(
CLOSE(12)
OPEN(13,FILE='VOR.DAT')
DO K=1,NN
WRITE(13,'(
计算结果见图2.2.3和图2.2.3：（这里给出的图的范围为20°N-50°N，100°E-130°E） 这里用的都是双精度，即REAL(8)，如只用单精度，将其改为REAL(4)即可。
2.2.7附注
图2.2.3 500hPa散度场（单位：10秒）
-5-1
图2.2.4 500hPa涡度场（单位：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值，单位：