第二章诊断分析2(2)

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