2.3.6例
以1988年5月1日60°E—180°—160°W，0°—70°N范围内的u、v场、温度场，用运动学法和求解准地转OMEGA方程法计算垂直速度ω，资料共7层，网格点为2.5°×2.5°。计算中采用插值法，将其插到10层等压面上计算。
PROGRAM MAIN
INTEGER,PARAMETER::L=57,M=29,N0=7,N=10 REAL(8),DIMENSION(L,M,N0)::U0,V0,T0 REAL(8),DIMENSION(N0)::P0
REAL(8),DIMENSION(L,M,N)::U,V,T REAL(8),DIMENSION(L,M,N)::OMEGA REAL(8),DIMENSION(N)::P REAL(8)::F0,DL,DM,PI
DATA P0/100,200,300,500,700,850,1000/
DATA P/1000,900,800,700,600,500,400,300,200,100/ L1=L
PI=3.1415926 F0=0*PI/180. DL=2.5*PI/180 DM=2.5*PI/180
WRITE(*,'(\，用运动学方法；IO=2，用OMEGA方程)\ READ(*,'(I1)')IO
OPEN(12,FILE='U880501.DAT')
READ(12,'(
OPEN(13,FILE='V880501.DAT')
READ(13,'(
OPEN(13,FILE='T880501.DAT')
READ(13,'(
CALL SPLINECALCULATE(P0,U0(I,J,1:N0),N0) DO K=1,N
CALL SPLINEEVALUATE(P(K),U(I,J,K)) END DO
CALL SPLINECALCULATE(P0,V0(I,J,1:N0),N0) DO K=1,N
CALL SPLINEEVALUATE(P(K),V(I,J,K)) END DO
CALL SPLINECALCULATE(P0,T0(I,J,1:N0),N0) DO K=1,N
1
CALL SPLINEEVALUATE(P(K),T(I,J,K)) END DO END DO END DO
IF(IO==1)THEN
CALL YUNDF(U,V,P,DL,DM,F0,L,M,N,OMEGA) OPEN(20,FILE='OMEGAY.DAT')
WRITE(20,'(
CALL OME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA) OPEN(20,FILE='OMEGA.DAT')
WRITE(20,'(
! SPLINE.FOR -- NATURAL CUBIC SPLINE
SUBROUTINE SPLINECALCULATE(INX,INA,NUMX ) IMPLICIT NONE
REAL(8)::INA(*),INX(*) INTEGER::NUMX
INTEGER::MAXPOINTS
PARAMETER(MAXPOINTS=1000 ) INTEGER::NXS
REAL(8)::X(0:MAXPOINTS)
REAL(8)::A(0:MAXPOINTS),B(0:MAXPOINTS) REAL(8)::C(0:MAXPOINTS),D(0:MAXPOINTS) COMMON /SPLINEDATA/ A, B, C, D, NXS, X INTEGER::I
REAL(8)::ALPHA(0:MAXPOINTS), L(0:MAXPOINTS) REAL(8)::MU(0:MAXPOINTS), Z(0:MAXPOINTS) REAL(8)::H(0:MAXPOINTS) NXS=NUMX-1 DO I=0, NXS
X(I)=INX(I+1) A(I)=INA(I+1) END DO
DO I=0, NXS-1
H(I)=X(I+1)-X(I) END DO
X(NXS+1)=40000 DO I=1, NXS-1
ALPHA(I)=3.0*(A(I+1)*H(I-1)-A(I)* &
2
(X(I+1)-X(I-1))+A(I-1)*H(I))/(H(I-1)*H(I)) END DO L(0)=1.0 MU(0)=0.0 Z(0)=0.0
DO I=1, NXS-1
L(I)=2.0*(X(I+1)-X(I-1))-H(I-1)*MU(I-1) MU(I)=H(I)/L(I)
Z(I)=(ALPHA(I)-H(I-1)*Z(I-1))/L(I) END DO L(NXS)=1.0 Z(NXS)=0.0 C(NXS)=0.0
DO I=NXS-1, 0, -1
C(I)=Z(I)-MU(I)*C(I+1)
B(I)=(A(I+1)-A(I))/H(I)-H(I)*(C(I+1)+2.0*C(I))/3.0 D(I)=(C(I+1)-C(I))/(3.0*H(I)) END DO END
SUBROUTINE SPLINEEVALUATE(EVALX,EVALY ) REAL(8)::EVALX, EVALY INTEGER:: MAXPOINTS
PARAMETER( MAXPOINTS=1000 ) INTEGER::NXS
REAL(8)::X(0:MAXPOINTS)
REAL(8)::A(0:MAXPOINTS), B(0:MAXPOINTS) REAL(8)::C(0:MAXPOINTS), D(0:MAXPOINTS) COMMON /SPLINEDATA/ A, B, C, D, NXS, X REAL(8)::TERM INTEGER::I I=0
DO WHILE(.NOT.((X(I).LE.EVALX).AND.(X(I+1).GT.EVALX))) I=I+1 END DO
TERM=EVALX-X(I)
EVALY=A(I)+B(I)*TERM+C(I)*TERM*TERM+D(I)*TERM*TERM*TERM END
计算结果见图2.3.1和图2.3.2：（这里给出的图的范围为20°N-50°N，100°E-130°E）
2．4通量计算程序
2.4.1功能
给出通量计算的程序。
2.4.2方法说明
3
2.4.3子程序语句
SUBROUTINE TONGLIANG(L,M,N,U,V,OMEGA,Q,QX,QY,QP)
L——输入参量，整型变量，东西方向格点数。 M——输入参量，整型变量，南北方向格点数。 N——输入参量，整型变量，铅垂方向层数。
U——输入参量，实型三维数组，大小为（L*M*N），风速的东西分量，单位：米/秒。 V——输入参量，实型三维数组，大小为（L*M*N），风速的南北分量，单位：米/秒。 OMEGA——输出参量，实型三维数组，大小为（L*M*N），P坐标下的垂直速度，单位：hPa/秒。
Q——输入参量，实型三维数组，大小为（L*M*N），比湿，单位：千克/千克。 QX——输出参量，实型三维数组，大小为（L*M*N），东西方向单位面积的水汽通量，单位：千克/（秒·米·百帕）。
QY——输出参量，实型三维数组，大小为（L*M*N），南北方向单位面积的水汽通量，单位：千克/（秒·米·百帕）
QP——输出参量，实型三维数组，大小为（L*M*N），垂直方向单位面积的水汽通量，单位：千克/（秒·米2）。
2.4.4哑元说明
2.4.5子程序
SUBROUTINE TONGLIANG(L,M,N,U,V,OMEGA,Q,QX,QY,QP) REAL(4),DIMENSION(L,M,N)::U,V,OMEGA,Q,QX,QY,QP REAL(4),PARAMETER::G=9.8 DO I=1,L DO J=1,M DO K=1,N QX(I,J,K)=Q(I,J,K)*U(I,J,K)/G QY(I,J,K)=Q(I,J,K)*V(I,J,K)/G QP(I,J,K)=Q(I,J,K)*OMEGA(I,J,K)/G END DO END DO END DO END
３．5流函数、速度势函数计算程序
2.5.1功能
给出通量计算的程序。
2.5.2方法说明
2.5.3子程序语句（包括流函数LHS和速度势函数SHS）
SUBROUTINE LHS(L,M,N,U,V,P,F0,DL,DM,FAI,KAP)
L——输入参量，整型变量，东西方向格点数。 M——输入参量，整型变量，南北方向格点数。 N——输入参量，整型变量，铅垂方向层数。
U——输入参量，实型三维数组，大小为（L*M*N），风速的东西分量，单位：米/秒。 V——输入参量，实型三维数组，大小为（L*M*N），风速的南北分量，单位：米/秒。
4
2.5.4哑元说明
P——输入参量，实型一维数组，大小为N，铅垂方向N层的气压值，单位：百帕。 F0——输入参量，实型变量，南边界的纬度。
DL——输入参量，实型变量，东西方向格点间距。 DM——输入参量，实型变量，南北方向格点间距。 FAI——输出参量，实型三维数组，大小为（L*M*N），流函数。 KAP——输出参量，实型三维数组，大小为（L*M*N），速度势函数。
2.5.5子程序（包括流函数LHS和速度势函数SHS）
!
求流函数和速度势函数，采用求坐标系，使用网格点资料。 SUBROUTINE LHS(U,V,P,L,M,N,F0,DL,DM,FAI,KAP)
REAL(8),DIMENSION(L,M,N)::U,V,VOR,DIV,FAI,KAP !,U1,V1,U2,V2,U3,V3 REAL(8),DIMENSION(N)::EMS,P
REAL(8)::ZM1,ZM2,ZM3,ZM4,ZM5,ZM6,ZM7,ZM8 REAL(8)::F0,DL,DM,ALF
REAL(8),PARAMETER::AA=6.371E6;PI=3.1415926 ALF=1.8323 L1=L-1 M1=M-1 DO K=1,N
ZM1=-0.5*(U(1,1,K)+U(1,M,K)) !西边界 ZM2=0.5*(ABS(U(1,1,K))+ABS(U(1,M,K))) DO J=2,M1
ZM1=ZM1-U(1,J,K)
ZM2=ZM2+ABS(U(1,J,K)) END DO
ZM1=ZM1*DM ZM2=ZM2*DM
ZM3=0.5*(V(1,M,K)+V(L,M,K)) ZM4=0.5*(ABS(V(1,M,K))+ABS(V(L,M,K))) DO I=2,L1
ZM3=ZM3+V(I,M,K)
ZM4=ZM4+ABS(V(I,M,K)) END DO
ZM3=ZM3*COS(F0+(M-1)*DM)*DL ZM4=ZM4*COS(F0+(M-1)*DM)*DL ZM5=0.5*(U(L,1,K)+U(L,M,K)) ZM6=0.5*(ABS(U(L,1,K))+ABS(U(L,M,K))) DO J=2,M1
ZM5=ZM5+U(L,J,K)
ZM6=ZM6+ABS(U(L,J,K)) END DO
ZM5=ZM5*(DM) ZM6=ZM6*(DM)
5
!北边界
!东边界