其次章诊断分析2

本文格式为Word版，下载可任意编辑——其次章诊断分析22.3.6例

以1988年5月1日60°E—180°—160°W，0°—70°N范围内的u、v场、温度场，用运动学法和求解准地转OMEGA方程法计算垂直速度ω，资料共7层，网格点为2.5°×2.5°。计算中采用插值法，将其插到10层等压面上计算。

PROGRAMMAIN

INTEGER,PARAMETER::L=57,M=29,N0=7,N=10REAL(8),DIMENSION(L,M,N0)::U0,V0,T0REAL(8),DIMENSION(N0)::P0

REAL(8),DIMENSION(L,M,N)::U,V,TREAL(8),DIMENSION(L,M,N)::OMEGAREAL(8),DIMENSION(N)::PREAL(8)::F0,DL,DM,PI

DATAP0/100,200,300,500,700,850,1000/

DATAP/1000,900,800,700,600,500,400,300,200,100/L1=L

PI=3.1415926F0=0*PI/180.DL=2.5*PI/180DM=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!垂直插值DOI=1,LDOJ=1,M

CALLSPLINECALCULATE(P0,U0(I,J,1:N0),N0)DOK=1,N

CALLSPLINEEVALUATE(P(K),U(I,J,K))ENDDO

CALLSPLINECALCULATE(P0,V0(I,J,1:N0),N0)DOK=1,N

CALLSPLINEEVALUATE(P(K),V(I,J,K))ENDDO

CALLSPLINECALCULATE(P0,T0(I,J,1:N0),N0)DOK=1,N

1

CALLSPLINEEVALUATE(P(K),T(I,J,K))ENDDOENDDOENDDO

IF(IO==1)THEN

CALLYUNDF(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

CALLOME(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)ENDIFEND

!SPLINE.FOR--NATURALCUBICSPLINE

SUBROUTINESPLINECALCULATE(INX,INA,NUMX)IMPLICITNONE

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,XINTEGER::I

REAL(8)::ALPHA(0:MAXPOINTS),L(0:MAXPOINTS)REAL(8)::MU(0:MAXPOINTS),Z(0:MAXPOINTS)REAL(8)::H(0:MAXPOINTS)NXS=NUMX-1DOI=0,NXS

X(I)=INX(I+1)A(I)=INA(I+1)ENDDO

DOI=0,NXS-1

H(I)=X(I+1)-X(I)ENDDO

X(NXS+1)=40000DOI=1,NXS-1

ALPHA(I)=3.0*(A(I+1)*H(I-1)-A(I)*PI=3.1415926ALF=1.8323L1=L-1M1=M-1DOK=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)))DOJ=2,M1

ZM1=ZM1-U(1,J,K)

ZM2=ZM2+ABS(U(1,J,K))ENDDO

ZM1=ZM1*DMZM2=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)))DOI=2,L1

ZM3=ZM3+V(I,M,K)

ZM4=ZM4+ABS(V(I,M,K))ENDDO

ZM3=ZM3*COS(F0+(M-1)*DM)*DLZM4=ZM4*COS(F0+(M-1)*DM)*DLZM5=0.5*(U(L,1,K)+U(L,M,K))ZM6=0.5*(ABS(U(L,1,K))+ABS(U(L,M,K)))DOJ=2,M1

ZM5=ZM5+U(L,J,K)

ZM6=ZM6+ABS(U(L,J,K))ENDDO

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)))DOI=2,L1

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

ZM8=ZM8+ABS(V(I,1,K))ENDDO

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

!南边界

EMS(K)=-(ZM1+ZM3+ZM5+ZM7)/(ZM2+ZM4+ZM6+ZM8)!订正系数DOI=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))!北边界ENDDODOJ=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))!东边界ENDDOENDDODOK=1,N

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

FAI(1,J+1,K)=FAI(1,J,K)-(U(1,J+1,K)+U(1,J,K))/2*AA*DMENDDODOI=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)ENDDO

DOJ=M,2,-1!东边界

FAI(L,J-1,K)=FAI(L,J,K)+(U(L,J-1,K)+U(L,J,K))/2*AA*(-DM)ENDDODOI=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)ENDDO

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

FAI(1,J,K)=FAI(1,J,K)-FFF*(J-1)ENDDODOI=2,L

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

DOJ=M1,1,-1

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

6

DOI=L1,1,-1FAI(I,1,K)=FAI(I,1,K)-FFF*(M-1+L-1+M-1+L-I)ENDDOENDDOCALLDIVER(U,V,P,DL,DM,F0,L,M,N,DIV)CALLVORTICITY(U,V,DL,DM,F0,L,M,N,VOR)VOR=-VORKAP=0

CALLERROR(FAI,VOR,L,M,N,F0,DL,DM,ALF)CALLERROR(KAP,DIV,L,M,N,F0,DL,DM,ALF)ENDSUBROUTINEERROR(A,Q,L,M,N,F0,DL,DM,ALF)REAL(8),DIMENSION(L,M,N)::A,QREAL(8),DIMENSION(L,M)::A1,A2,RREAL(8)::F0,DL,DM,AA,ALF,EPSAA=6.371E6EPS=1.E2L1=L-1M1=M-1DOK=1,N

WRITE(*,'(\N0=0DOA1(1:L,1:M)=A(1:L,1:M,K)A2=A1N0=N0+1DOI=2,L1DOJ=2,M1R(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)ENDDOENDDOIF(MAXVAL(ABS(A1-A2))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=