其次章診斷分析2

本文格式為Word版，下載可任意編輯——其次章診斷分析22.3.6例

以1988年5月1日60°E—180°—160°W，0°—70°N范圍內(nèi)的u、v場(chǎng)、溫度場(chǎng)，用運(yùn)動(dòng)學(xué)法和求解準(zhǔn)地轉(zhuǎn)OMEGA方程法計(jì)算垂直速度ω，資料共7層，網(wǎng)格點(diǎn)為2.5°×2.5°。計(jì)算中采用插值法，將其插到10層等壓面上計(jì)算。

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(*,'(\，用運(yùn)動(dòng)學(xué)方法；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