第二章诊断分析2

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,'(F7.1)')(((U0(I,J,K),I=1,L),J=1,M),K=N0,1,-1) CLOSE(12)

OPEN(13,FILE='V880501.DAT')

READ(13,'(F7.1)')(((V0(I,J,K),I=1,L),J=1,M),K=N0,1,-1) CLOSE(13)

OPEN(13,FILE='T880501.DAT')

READ(13,'(F7.1)')(((T0(I,J,K),I=1,L),J=1,M),K=N0,1,-1) CLOSE(13) T0=T0+273.15 ! 垂直插值 DO I=1,L DO J=1,M

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,'(E11.4)')(((OMEGA(I,J,K),I=1,L),J=1,M),K=1,N) CLOSE(20) ELSE

CALL OME(U,V,T,P,L,M,N,F0,DL,DM,OMEGA) OPEN(20,FILE='OMEGA.DAT')

WRITE(20,'(E11.4)')(((OMEGA(I,J,K),I=1,L),J=1,M),K=1,N) CLOSE(20) ENDIF END

! 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

!北边界

!东边界

ZM7=-0.5*(V(1,1,K)+V(L,1,K)) ZM8=0.5*(ABS(V(1,1,K))+ABS(V(L,1,K))) DO I=2,L1

ZM7=ZM7-V(I,1,K)

ZM8=ZM8+ABS(V(I,1,K)) END DO

ZM7=ZM7*COS(F0)*(DL) ZM8=ZM8*COS(F0)*(DL)

!南边界

EMS(K)=-(ZM1+ZM3+ZM5+ZM7)/(ZM2+ZM4+ZM6+ZM8) !订正系数 DO I=1,L !对边界上的U、V进行订正 V(I,1,K)=V(I,1,K)-EMS(K)*ABS(V(I,1,K)) !南边界 V(I,M,K)=V(I,M,K)+EMS(K)*ABS(V(I,M,K)) !北边界 END DO DO J=1,M

U(1,J,K)=U(1,J,K)-EMS(K)*ABS(U(1,J,K)) !西边界 U(L,J,K)=U(L,J,K)+EMS(K)*ABS(U(L,J,K)) !东边界 END DO END DO DO K=1,N

FAI(1,1,K)=0 DO J=1,M1 !西边界

FAI(1,J+1,K)=FAI(1,J,K)-(U(1,J+1,K)+U(1,J,K))/2*AA*DM END DO DO I=1,L1 !北边界

FAI(I+1,M,K)=FAI(I,M,K)+(V(I+1,M,K)+V(I,M,K))/2*AA*DL*COS(F0+(M-1)*DM) END DO

DO J=M,2,-1 !东边界

FAI(L,J-1,K)=FAI(L,J,K)+(U(L,J-1,K)+U(L,J,K))/2*AA*(-DM) END DO DO I=L,2,-1 !南边界

FAI(I-1,1,K)=FAI(I,1,K)-(V(I-1,1,K)+V(I,1,K))/2*AA*(-DL)*COS(F0) END DO

FFF=FAI(1,1,K)/(2*(L+M-2)) DO J=2,M

FAI(1,J,K)=FAI(1,J,K)-FFF*(J-1) END DO DO I=2,L

FAI(I,M,K)=FAI(I,M,K)-FFF*(M-1+I-1) END DO

DO J=M1,1,-1

FAI(L,J,K)=FAI(L,J,K)-FFF*(M-1+L-1+M-J) END DO

6

DO I=L1,1,-1 FAI(I,1,K)=FAI(I,1,K)-FFF*(M-1+L-1+M-1+L-I) END DO END DO CALL DIVER(U,V,P,DL,DM,F0,L,M,N,DIV) CALL VORTICITY(U,V,DL,DM,F0,L,M,N,VOR) VOR=-VOR KAP=0

CALL ERROR(FAI,VOR,L,M,N,F0,DL,DM,ALF) CALL ERROR(KAP,DIV,L,M,N,F0,DL,DM,ALF) END SUBROUTINE ERROR(A,Q,L,M,N,F0,DL,DM,ALF) REAL(8),DIMENSION(L,M,N)::A,Q REAL(8),DIMENSION(L,M)::A1,A2,R REAL(8)::F0,DL,DM,AA,ALF,EPS AA=6.371E6 EPS=1.E2 L1=L-1 M1=M-1 DO K=1,N

WRITE(*,'(\ N0=0 DO A1(1:L,1:M)=A(1:L,1:M,K) A2=A1 N0=N0+1 DO I=2,L1 DO J=2,M1 R(I,J)=((A1(I+1,J)+A1(I-1,J))/(DL*COS(F0+(J-1)*DM))**2+(A1(I,J+1) & +A1(I,J-1))/DM**2-TAN(F0+(J-1)*DM)*(A1(I,J+1)-A1(I,J-1))/ & (2*DM)+AA*AA*Q(I,J,K))/(2/(DL*COS(F0+(J-1)*DM))**2+2/DM**2) A1(I,J)=(1-ALF)*A1(I,J)+ALF*R(I,J) END DO END DO IF(MAXVAL(ABS(A1-A2))

2.5.6例

以1988年5月1日60°E—180°—160°W，0°—70°N范围内的u、v场，计算流函数和速度势函数，资料共7层，网格点为2.5°×2.5°。

7

PROGRAM MAIN

INTEGER,PARAMETER::L=57 !计算的东西方向的格点数 INTEGER,PARAMETER::M=29 !计算的南北方向的格点数 INTEGER,PARAMETER::N=7 !垂直方向层数 REAL(8),DIMENSION(L,M,N)::U,V,KAP,FAI REAL(8),DIMENSION(N)::P REAL(8)::F0

REAL(8),PARAMETER::AA=6.371E6 !地球的半径 REAL(8),PARAMETER::PI=3.1415926

REAL(8),PARAMETER::DL=PI*2.5/180 !资料的网格距(由经纬度换为弧度单位) REAL(8),PARAMETER::DM=PI*2.5/180 !资料的网格距(由经纬度换为弧度单位) DATA P/1000,850,700,500,300,200,100/ F0=PI*0./180 !南边界的纬度 L0=L

OPEN(12,FILE='..\\PRO\\U880501.DAT')

READ(12,'(F7.1)')(((U(I,J,K),I=1,L),J=1,M),K=1,N) CLOSE(12)

OPEN(13,FILE='..\\PRO\\V880501.DAT')

READ(13,'(F7.1)')(((V(I,J,K),I=1,L),J=1,M),K=1,N) CLOSE(13)

CALL LHS(U,V,P,L,M,N,F0,DL,DM,FAI,KAP) OPEN(15,FILE='FAI.DAT')

WRITE(15,'(E11.4)')(((FAI(I,J,K),I=1,L),J=1,M),K=1,N) CLOSE(15)

OPEN(16,FILE='KAP.DAT')

WRITE(16,'(E11.4)')(((KAP(I,J,K),I=1,L),J=1,M),K=1,N) CLOSE(16) END

8

本文来源：https://www.bwwdw.com/article/loxp.html